{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 12.8 Additional material\n",
    "\n",
    "This section contains additional material concerning confidence intervals for a fitted value and reference ranges. These are useful topics in regression, but will not be examinable. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12.8.1 Confidence intervals for a fitted value \n",
    "\n",
    "So far we have only discussed conducting inference on the estimated regression coefficients. However, it may also be of interest to determine **confidence intervals for the fitted outcomes**, or **prediction intervals**. The subsequent two sections describe and illustrate these two concepts, respectively. \n",
    "\n",
    "Rather than focusing on associations between variables and the outcome, we are sometimes interested in the expected value of the outcome at particular values of $X$, i.e. $y_x = E[Y | X=x] = \\beta_0+\\beta_1 x$. \n",
    "\n",
    "The fitted value is our estimate of this quantity, \n",
    "\n",
    "$$\n",
    "\\hat{y}_x = \\hat{\\beta}_0 + \\hat{\\beta}_1 x.\n",
    "$$\n",
    "\n",
    "The variance of the fitted value is given by:\n",
    "\n",
    "$$\n",
    "V(\\hat{y}_x) = \\sigma^2 \\left(\\frac{1}{n}+\\frac{(x-\\bar{x})^2}{SS_{xx}}\\right)\n",
    "$$\n",
    "\n",
    "where $SS_{xx}=\\sum_{i=1}^n(x_i-\\bar{x})^2$, i.e. the sum of squares of $X$.  \n",
    "\n",
    "The 95% confidence interval for the fitted value is given by:\n",
    "\n",
    "$$\n",
    "\\hat{y_x} \\pm t_{n-2, 0.975}\\hat{\\sigma} \\sqrt{\\frac{1}{n}+ \\frac{(x-\\bar{x})^2}{SS_{xx}}}\n",
    "$$\n",
    "\n",
    "* 95% confidence intervals can be obtained for values of the independent variable that do not arise in the data. However, the width of the confidence interval increases with the distance from the mean (as can be seem from the formula and figure given below). \n",
    "* Care must be taken when extrapolating outside the range of the observed data as this makes an un-testable assumption that linearity continues outside the observed data range. \n",
    "\n",
    "*Example*. The R code below calculates a 95% confidence interval for the fitted value of birthweight of a baby born after 280 gestational days. \n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Example\n",
    "\n",
    "We return to our first example, exploring the association between birthweight and length of pregnancy (gestational days).Suppose we are interested in the expected birthweight for a baby who is born at 280 days' gestation. \n",
    "\n",
    "The code below refits our model and uses it to estimate the expected birthweight for this gestational age. It also provides a 95% confidence interval around that estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>fit</th><th scope=col>lwr</th><th scope=col>upr</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>119.8818</td><td>118.9215</td><td>120.8421</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lll}\n",
       " fit & lwr & upr\\\\\n",
       "\\hline\n",
       "\t 119.8818 & 118.9215 & 120.8421\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| fit | lwr | upr |\n",
       "|---|---|---|\n",
       "| 119.8818 | 118.9215 | 120.8421 |\n",
       "\n"
      ],
      "text/plain": [
       "  fit      lwr      upr     \n",
       "1 119.8818 118.9215 120.8421"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Refit model\n",
    "data<- read.csv('https://www.inferentialthinking.com/data/baby.csv')\n",
    "model1<-lm(Birth.Weight~Gestational.Days, data=data)\n",
    "\n",
    "# Confidence interval for a fitted value \n",
    "new.data<-data.frame(Gestational.Days=280)\n",
    "predict(model1, newdata=new.data, interval=\"confidence\", level=0.95)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We estimate that the expected (average) birthweight for babies born at 280 days' gestation is 119.9oz. The 95% confidence interval for this estimate is is (118.9, 120.8). Informally, we can interpret this as: it is plausible that the true value of the expected birthweight, for babies born at 280 days' gestation, lies between 118.9oz and 120.8oz.  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can extend this idea to graph the fitted values - estimated expected birthweight - and their confidence intervals across the range of gestational days. The code below does this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAJYCAMAAACJuGjuAAAAM1BMVEUAAAAAAP9NTU1oaGh8\nfHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////UNI3wAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAgAElEQVR4nO2dibqjIAxGGbu3ty3v/7RTFUjCoqigoPm/mS6KxOXcJESr\nQrJYGSS2XgHWPsVgsbKIwWJlEYPFyiIGi5VFDBYrixgsVhYxWKwsYrBYWcRgsbKIwWJlEYPF\nyiIGi5VFDBYrixgsVhYxWKwsYrBYWcRgsbKIwWJlEYPFyiIGi5VFDBYrixgsVhYxWKwsYrBY\nWcRgsbKIwWJlEYPFyiIGi5VFDBYrixgsVhYxWKwsYrBYWcRgsbKIwWJlEYPFyiIGi5VFDBYr\nixgsVhYxWKwsYrBYWVQpWO+rEKdH//krlNovd9Hc1dTGXujvehKiuTyjLNwbIa6//SPIHrK+\nzlKoj2vS3rZWmWs1pldP0rn78ofA6ma82qk38bAWuuh2zd+4hXvbcE2w/pp5PTNYCfXWiFza\nb3cE1kX8/fVTG/GlC50F6D1q4qQarQbW3J4ZrIT6xcHLV77Vwf8RY1xQu5e7Pf0QN7rMz181\njx9rn4cGclCB45XvMDJYBeiX/rTu6Oe42nwK71oD1kl8yCK/to2a8j7d+0+vaxvuXnpJIV8/\nRq9v8xW/f65Nm9Tpr99bI5rbx7fkb96P+PNLOg3NOgpnMZQlWl1/TuL210Vl2f1FtX9Dzzaq\nn/o2eql765Ej88c1VCVYcNQvfVy86sOqwXrZTukmnJxLx8aL7vNm4iQcZ2VK5XFn9fXT9N//\nnCVh3s1pSNc+YNDuuoW0/UtSizZoxbs2dCmVdhagWsFqPda326kPvZtbH6FzrLOwEvTfwaAu\nDHL5nixIwK4esBqYjb82zpKo6ctuCGs/YNDT9bP9s2h7e3W4/jb4/O3+VC6mt2vbSn7P7p/P\nVqoSrGu3a7+9/wBAPmZU+BYnaxFhpyJty1/O9b0rBNpj+eq6Frh9//7sZ74aoRK4zn5/NK0l\nf/Oad7duJ7shXRW/Qafrswr6V6kj4Un9jWDy+61v/9TsDd9KVYJlRoXtTr2eW0De5z743Js2\n77qKl7WIA9ZV/3Hf+oOm+Pp6wbqoGsbLfP32My/Okqrpt0vkaEO6Kn6D3q47mL42Nhislvir\nvc2bqkqwdB2rQbS80V7viqP3RpVKOzlgCXUE5YfEPOkFyyxtvqKAFWgqnYZ0VfwGna77lby3\nHuzZD1badX7ezgLbUyWXgtiqE6yu8n55nzzHq1VbHIVSaaeTXbyC5kvAcpcMgEUnuhZcsEhn\nnzYt13ni8+Q2uWkcrUxyM1UKVq8uXlwvtDAg++IolEo7OaNC7LFcvzMMFvaUg2DRhrS936C/\n627woX3yz3GJ0/Xxpva+z36wWMqwsEqwLpcOilcHy7l3TH+wT7viqKlo9fqDP+Y/nf7YOVb3\n1XvYdY71pF+ldJuerRzLDU6DYPm77izrNT6pJi7Ir6vwkbyJSlmPSbp0LLz6Omk/DHs3MPLq\niqMWWFB5b9ORi3dU2LXzHvZHP4J7qlHhsz/f+OxZdpuaUSFtqBQG6xvqWp9p/6IFiMc6mXnO\nqfeNVCVYL51jdLmsrvyczNx+QEXB+qJSVOe74NxhV34aBCtUx0IlSrfpw26oFLZwC3Ut+7JE\nT2c//tWlj/71t0fOny6Ht85kbaYqwdK1qz6F+kO0tDp1x8TKsfrqoSKwb6m/92dLhsFSNi7q\nqyb75muKK++koZLfggbH27Wa/sTb246J/5zkvZQUq1Kw5OPn+y86GWnP4zU3fTGDKo5ao8Ju\n0q/dL+19ke/4XCF6t77+bPwOGj5XCCtgN21P9pmVww0lsWQvdlGI+7pWn9UmtkPi5vr+dH9Z\nukmXX51LqbtXC9aALgqnvlTK2kj7A4tVhBgsVhYxWKwsYrBYWcRgsbKIwWJlEYPFyiIGi5VF\nDBYrixgsVhYxWKwsYrBYWcRgsbKIwWJlEYPFyiIGi5VFDBYrixgsVhYxWKwsYrBYWcRgsbKI\nwWJlEYPFyiIGi5VFDBYri+aD9Xfv76BwuUU86IF1NM0F62vuKlfSnShYxWguWDfRPPt7L35e\nTTH3zmEVo7lgNeiWnu9i7vbFKkZzwbJulZliVVh7EnssVhYtyLFe/Y3xOMdieTQ7iOHH/52+\n4+1Zx9KCOtatq2M1lzvXsViOOO1mZRGDxcqi2WB9rt3z4h8n0XDqznI0+5ROdzfzx51P6bC8\nml9u+PmpWyOuX/m9DZcbBKt6rQdWo26C3z/BZbBAymlc9VoRLP10BfSFzl6EO6ssbeCx2tcv\ne6yda0WwdI51+6rP6U2witGKYE0YFTJY1WtFsCbUsRis6rUmWEWZYOUVg8XKIgaLlUUM1o5U\nUgGQwdqN1IPEt14NpVUr79HF9VL2TlUS6HV7rQjWg8HKKWG9b6w1Q+G7ib1YppCdU5UODJZ8\nx/42p5CdU5WODNYvGr7HGy0zcQQJlacL9O3AOVZZJuoVour3AmPBw44KyzJRr6h/It+KwYrB\nqk8+otQ7g1WAiWoVBItDYREmqlUYLDJhazFY1SmQYx263FCQiXoVGBUyWGWYqFneOhaDVYaJ\n/YlzrCJM7E88KizCxB7lq2NtVdtisHat7bwYg7VrbZd3MVi7Ewp+G44UGaydiQQ/Bmt9E3sV\nCX4M1vomdioLJc6xVjexU9lg8ahwbRM7lRP8uI61son9iKJTyokdBqty2cGulBM7DFblcj1U\nGRcoM1h1K1hQ2BovBqtuBcDaPiAyWHUrBJZv4qpisKoUBDovQgVcTcpgVSgc6LxBj8HazkTN\nol7Kk6YfDay/e/eEVXG5jTxhlcEaEsImNPY7VI71PaHbrvEDBOaLuiovWocaFd5E8+zvYjT6\nFHsGa0gAFvlqNTpOHatBN8d680OaZkvgtF1NsHfY1litChY9VzrYy9Z7pWAJnVkZdtywt30g\nZI9VnbS3EvSaK3q2UNpT1te6Odbr033iHGu23BGhfX17EcWGdcsNZzQqPH2zmNi9yIhQvdA5\nBwRL/t26OlZzuXMda6bsEaF9boe02TKJ58p7ZbJHhL6MSk/ZMolnsCoTuYOR8kpS+keFWybx\na4L1vbVDwftJiPMzk4ndyXci0LnNmuhuwIZoE2oKarO2VgTr0/y2tH8wNJ/SidNALHNORAud\ncBmfJjyViNW0IlhXcfn+Xq6f7unQXG6I0EAso8z1YU+/qoW0Y9s7WEJ81csvKnKBNELDdQNy\n8w+Tryu/JYCs3edY3X5oSLmYzo585txxFF2QQmChQChROFxdq4bC9y9178/rfIeTLAar02yw\nzJDwEHWst2hub3lpfmS9TuKVw8R+pDxP/3m4kXFQEo0NBy+pWUNrlhteDcS6ex4TO5H7aK+R\nRmZUCCd6NiyOdmu3yiJaz2t3Fenl/slmYhcCdzMQyyBflxJqV/rzlmGwX681FinQRMmKyq2s\nRrpsigPhivpnT2CwCtR0sNT4z3uCZwX9c7hisErUDLCkOWvou1I5s1ysGKwyFTWkQzmWICNC\nyLryrSGWh6tZthms7IqKZ2hUiAujeOk10PKEwZn1/vmV9+ji+tHBivyTt+tYEl3TsFYdy+uu\n5lmeu7IPBiuBMEzShD9dT3Uulcm9J/3uap7h2ev6boYvlklgYu+yKqemLGrKo/alMnn3pG80\naJE9QfPX9T18sUwKEzuXVRSFi2WEPu8Mw0SJ37Mo6K7mGV6wrg/008JMJnYmmjNQXAxRQo8K\nEVQr5FhhdzXPMo8KV5M9tHPAgitH7TCYfVQ4UhRlsEqW7XZcsFAYpB4rdx1r0F1JBqtouYmS\nlWNBvi6sUmnuVRs9h8NgFSwPWIFRoaILl0pzrpgbBh17DFbB8g3tSFHUZO5OqTTrLow55cxg\nlaxhRqCAhQsQqJCUJyRGuCvJYJWtkStI8XkM99eEmUKihyvv2k3vmcFaUUNXkOKyuzRh0JRI\n84TEKHc1zy6DVYRMDLR/dqHgylJ7j3RX88wyWEUoDBZMNdMSyeFq7Ir8SWKwVpfvchByoYgp\nwOtZ/iHlQk25sJ3BKl+CYgMTafIO8c99TaEJ7mqeXQZrZeE7fuCJFCxpBoGmpiVHDv4kuVyN\nrPRkMVjryrpIxhTdJeUKnBgsmbDYYHM1euH0dBMM1rqit5KRUK3C555Fcg9F9M/hamwJBqt4\nEbDUJw2UiYJ5wZrqriSDVYMgxxIm3gFZNERm2XcR7op/CV2hYFQIZ5xFQBn2nR0GPe6Kfwld\nqXRODok7HhVK8GMywzmcUXfl+2Uhg1WTtEfCLAmThNFSvGm/TOPuyosVg1WRIFn3XOYnyakc\n62W+HK7GGoRbjovB2kbaWwnsraS+MBk8lmqbovI+113Ns8tgZZSOXm4UA3jI9ckQH/Vc1Zqe\nKwxHxfCcCHc1sCHhWQkXKdBEkbIuaafz+ga42m7KpRAG4SdgGKxwVByIl2PuKhgFYXWnicHK\nJuptfGCZ6hVUtQxy2mnhPsaiYnjOEncV6HJEDFZimWBEo5eHLCgpGKJQCV5XG4RFkz+8wiKO\nJem4I4+7GtuokflpFinQRDFCwWgELFMntQql9hwcSoWZgOIk7cy15HI1PNu7WWMNkixSoIli\nhILRMFjaYYEPQiND/I3WsUxeb1yZxIv4LC12V551j9D8o/53756wKi634z5h1R6EIYjsER3h\nQ4cuARm8GSSihF7aO4/4NzPXmPAt8uMKr+QcrFYF63tCe+Cgj5VzB2HUu+gRHmpmvJAYFgBo\nWSRR0wsWXaT3V8aXWcciIgqS7ZqguUf9Jppnfxej4z7F3h2EAVjqzU60UQ6O3JVEl2QhsGx/\naPXrgiXtRVp3FVhVGY/VqmA16OZY72M+Vs6XPCG/ZM+SGBiVKLlhEDkkRInADg8CIWAGyRmW\n4UqzbM2dvKkTNPeoW4PdHCZKlxcsK+F25+Ahn3OpuzUJR1vUs0bOvKBesVqPRPN7a+bkTZ0g\n9liz5XVLZvBmz9I+BqfedtyzKVOLQRwzEDn5ux8rJzjTmZM3Ne8inX451qt/OBPnWDGz0PEX\nBCOcykuElOYzXHUX2ISzFhqdQHY1YTND25h8kV5ntI9O3ywmSpdv3Baa5YI1ODQMgoUDIvTr\ny8v/wZpY3myiu3J7z7SI0t+tq2M1lzvXscZnYbBIDAwWsqQPLJRm4Yn2TiboLMZqZbBKMlGB\nIKsarmGZxpgt00GgX1pA+2dxRdrP4YrBKljWqDAYElVj4rjCIZfO7D9TdBJgtRFYw7WGJCZK\n10BEpK1Imk5DH6qI4jYRBlCz9mXAXcWcbvZaYrC20KBDsdtKQ5XJlcyFpOiVOKvoFZF9GBRk\nCmicq5DNFcHyZQaJTdSiwRTIbmvtNAyWxgvO+EwHq2PHX7uKiYKhTVkRrL+GwerlHZcFd4md\nX9lZluYKUzbcHfryjxRFbaxiufIcsTVD4fcizl2F1B+UY6mrX56jEY5ipqqEaLL3Ft11YqA7\nawbhak6NoQiwpHwK8ZScY/nAsifAHIh1XpSkMFeHIrBC3ZEZOruihVM9a+6mBCbE9zVHn7O4\nfA8JFo1A6FVn4GQSaUrOHIdkl+j93enpQmftZsXIMYmKgr5NcSxNEVnk79aeqDmPXRIKuovm\ndTywrAjkVpJCJFhgDcNFIPN2pxlqmxB6Zrkrz5b5O4zsCj4+4aLQ0yty8fdJHBAs9Np/pE5q\nECxTr3Lg8qdaaIToXY8uVFKsZrora1McS9NkFvnFtfPj3Z5N/v7dzyoxj9D1cGCFIx2ctwnM\nN+dqMD6h5F1VtKSuRkjnsAuLKyEENTsVq5AWgPUSN3yJwucmYp1WtIm9aBgsNO7zLTtVpqgl\nPD+p7mf8U2m7M/tfKq6WgHWxr3z5XpeujW1iLxoBq3sPJwjgisJuCmOFPJYntRY6u1KW8YWi\nybBafVRYjol15T/A6l3NstAy3yDwGWB8RKHzPvq/NK4QG/6Huep9FsxIvMW5FynQxLryhSTr\nhTZB3ybHQkqhnWZproSkHi1dFIQ+5y9ifiAYUSx3NjzpWpUua4vRETWHH6Zb80nRPYQRbUAS\nd2S6x0fQ2pWQid2VTACWIisCrMehwaKitYbuRWCy0HwFC64x+JItqTky/soHFnFLyGEldlcy\nAVjXnqwIsOS7Gf7986K1qksULDwqtMAKuiiPzzJ9SISfxM4S84MatxPTYpUALHkWVxkHlnwP\n/zZn0VrVJWuYiNyVDZaUZEhoMifqsyD+QcrmJO+Wu9KW/6V3VzIFWD+ybpFg/aLhe7zRvLWq\nTNYwUR98TYQZ/vUzcQAMe7G+dAVDSKvcYLsrCVOTY5UErI6sOLDmmNinrGEijmFAlfZjlJ6h\ncGiFRmSGAESwSh4FqYl5i/Rb3ogbgzVVAu8yAEh91Y7JtB2kikKIfJYxs667konA+jSCwVom\nMnyDH8RDKm7C26j0oojdtd2VXAyWUkvW8nUZNLFzkeEbSbSgMjUQ+WgTnGbZYZCMECf/cH7C\nBq2ySIEmipPAxQYT10zpQFrDw7CvAqfV9yWD7iofVgnAerYX+l2eiVbHa6JwxaYBUe1wYqQv\n2NORbaT4jmoRUL+g1QTLXaXYquENmb+IvtFHbOlzhomyZQ3vlrZzRoVCVxqkHHNY9tyus3/o\nSlFSJx2+/X/kVg30sGyRh2jai7BejXgsWIlBE4XLKkgtbmfin9DF0UDOFXBZqJAqxVB2lWZt\nR3qYvchJVTzf4rRgJQZNlC1hvc9th8iRCC76MVx2IF5OaC8n/6GfOWs6Y5L22K2K6GLmIuav\n4KijwjRguWXR8ag3JnWdqLr3oyb2X1SNYXuwwGMN3vpxiYmyNQksdXRDvZixn1SNHbxkBGwm\nYP6DnzmbAWdkjWF7sDjHis+xTGAL9oHeNVXYg0W5KS1I0DWvMUn71K0a62H+IjwqjB4V6oTa\nM0+9hsAyX+PB+qd/L6FZVol8ZKV981GhlM/27o9cxxqfJswRdlurV3hHfsomapQwMxj8B6MB\n463+xe5a31ZN0WKwsqgqsFx5/96H8hY7x1qYvmuGhPGTxl2ttmcZrAzyZiiDYFmjQhQII9lC\nzYxrMoFXT/qXePA+pAVg8e8KA7IQ0kFlMCEWLlv0c7TzAq4kiYK/UeJ6WC0Ci38JHRABC8Ji\n/EmdaSR5sVJdtB0ad7WqloTCufduyLJWBYmChSZFuQw9eLQ8lXS+WJfBe7iSYiusluZYc+42\nM9FEjUIwWYVRCy2bNNVQMeFxSJ6JLlb/YBywFVbLk/fJ98eabqI+oaAnBoKhHRsBKnNFg0XO\nsMBdSXT+ZuWk3WzMKosUaCKvwEWRVEtPkugDAks31Gw5LikCK/10VDpl/T3KYGUWOC6SewmD\nHp6snQ0uvkfiZQqguuU/zNX6u5TByiyUNAFBJAiaybriYF7iT+MAVhZn/2x81xKDlVvaNxGw\n+hn6q25o2pMz0AYfD2eWdzJJGU7jqZW1xGBlF6RS+pMZIpJxoYDfEmrnNeaypC8KoqRdCOQG\n1xWDlV3uqBBqD8I4NBQKgaqYWEgpsiboUhdYWW2rV1mkQBNrCo6qQKNDKFvpb3os2H91a1bO\nlDbo2cV2ApqUwrKy2javsojS3717wqq4jJW9dgaWI51N+b7iqb7A5y2JWu5KQhQklY8VtRAs\nsxea8UuTvye0k4YvDNw9WDTxMVm7tLzbmCx35QGtdrA+ESH8Jppnf4X8cZ9ib6TClPqiX+nU\nsRQLlRTwBGGwgj7BymobOH+RF9nK8Z9/NejmWCM/vtg/WL0ETbnUq8nhNR9x7gp5K31SSI8B\nVL/k9BG4SbF85OguvMRj4dB2Gj9ZKHyubsTErmWgMd9lfFn0n58r0d+6VveiR4WogkEt4Yrs\n7DOKvoVT5VgRYo9lyfUTlh/xeC39PuyudKlfUytMjYO+CvQK02dtiL3wiqPCX4716i/Z4hyr\nlbDecTaEGbI/e9yVRNmVOR9kToBbfVPLmq0FZDkbMrOn2Uf9jHbNyb6uOY2JmmQdDxRPegcz\nEAbD7oqEQThrBF9qAethEq2IJf9uXR2rudwPXsfqZIMFrypyBTyW7a7+2bUr7bDQ2SJwX5WA\ndUfbnlCHAIseS8jgjbdxJQfdldRuSedYxBCOjeXnWIl/Wu8zsV+RwRQNiOEo6HdX4h/k62a4\nCYZUz1B62NWoUH6vQpxfMQseAyw7WvVTpEQlB5sxg5UKe2QwCESaxa2+hWNZlFjHku1IbzAJ\nJ/o23cZe1JrEmjiKtFtRb5PclTCOCoGFSmSLAt2CzVmwyOUc/TOKWxs2v48m4tk7hwQLxRM/\nWP+Elyt0fZ+uuQucPll9r7U58xexNnt0uaZv8mlOHwbLFvUwpGw+7K7U5TGOt5Lm3fS/7gbN\nX2QqWLrJ93z2gTWts30JPIoJWg5X5KYMEAX/CVcySd60cJNWWaTTyeRjpzN7LCLIgQxiVhQc\nclcOU1ANPQRYD6FvGvIRZwYLiYzatMcmXEkCkuOu0CkgGBmastaW2zR3EbTx59FnEcKTnF4j\nf0pHBMt4G6gxTXVXyG8BaFtu0/xFyOaMXkT6vuhPnyuDBSKVcZRleUvt2F1ZObup28OHWsGS\nV31z2z95iX1+6kQTuxZyMMAXVA2QexpyV9ANLKh7U1bi1mVC67HOli1yM7fjPstvuocIHAUs\nCgXk24BMKAraYVB6vuqxd1wVC5/mSZCbpTqlA5uRQocBCy4uMEVynXW3Uy2OvFjpciieImEM\nIGPr7rpVoir9QrAa/ACBfGAlcc7lCV04Bd5FRzNJzgSG3RVNdJHnc66bGFmXKa1jtm3RIjfz\nAIGbfKa72Ttdqw1OSKwjYaqYwmylMPXNKKxwUmXBVjVY+AECYuQSGmfjY9cqkXMuTwQsSIt6\nVOK4smv00jitysGSr+4BAq3bEvfh5R4zwUq1qeWJpEMIj9go6PopPRww6RuUHOxdbh+EsnKs\nSXo3sbHyMGAJcDDms5DhGsNwTRTgBMCge5JOIMtmiixoVDhN79hK11HAksRLmbhoYxTESrrf\nTKcS5RAaIISRGY+i1SmkjmUNSGIWfaCfFsav1R5zLHTkdSlT1xkAI6tE6joon+OSulqB/yIF\ntiqBKJFtt64L1ry12t+o0BN2jMOKdlcB4TGBsgYEC2KtULAyyjaRmNvNBT7Y/EkqsDBGPsx8\nXgpPR3UGFyz8kcHakQAiid+h5uDcpDbaXaGqhUmqFLTKWwnElz/HSrutSxdpyw1SXtI97sQ1\nsRNB3LPB0h7HwWhqFEQJljPkhJo8sZltaxcuclZr3iQla59gmVdnnNsfYEjaB7J2hyWMlQZL\ndWp8JMm5EM4Z9/NCsB7i/G1XD64OTaI9goUzHivF7mYhjChX1OcEudJfjWsEqAxKeD1W2tyZ\nizTim2PQtmuwoHyJ2JK2d/K7q1DuHvZoNlhrDYMWgoX+QFKtkW1iL8JVGuOv9BRNkdQVHAcr\ngpSkIZD4K4kcnCRgZc2pQps7c5GT8ljvdBf52SZ2I1wE0AdcTdM/llANBPJfUgoPQ/TEovFN\nqHoPcEmd0K+6V9PkWK/ENwfZJ1gohSZgYYq6BlYty0OWL92iYJlWyPSqW7twkYvarGSXYrkm\nEnVZQI1V5zwGrJ4iicOgdItZUdKeTSr/pBP43vCIx8qwcxaC9acum3kmWyHHRJoO1/+TDcr4\nkO5w0uxKwrlCh5yAp1J1BvwNUjiS0IW2P8vOWZq8N/e0pVHXRMIOCwELZdgIKzWDuis7Rw9E\nRjpF4uqE6TjMTpadsxCs62/Vz8/4WxnNMJGyvw3Jwv5CRy76yGacb3mjnUMTjAclnqBMIMy0\nXf+aDc6dvb1LF3m2tfdryieNOybS9bcZWL5oY7sr5K8GmPL5LDWU1BUGA5OZNeixCgVLys/9\nJEST7MeqPhOp+tsOLNe8wQqyeQNalL9C00zJypSu0AkdFB4H1q1AsPp7QCZdrb3lWJ5jhyFS\nvgQ5sAGybM40TMCXGR8IbVU3HVq7snKsn96twxLnkV9SLDKRoMNtR4UOWAQrkrXjFGnUVeF6\nFbgvGxXTJLR2BY4KX7dGiNMtcYq1uzqWBZaDlbC4siEajI0mqwJ2hGN1ePMz7Jyl5QYhLnGX\nsc82sQuRCiUOeRouVM3qpnsinjQnagA16APT4ydq3Z261GO12dXPYyUuOOwPLIg2eCio02qn\nmuX3U56yldQDRVJ6MFZxDKwKLNk+x6Tp4EqzPl4T9ct4LIKVl7Vg5g4FBQKaRCNAIa2opp1d\nvw5rKsmo8K/4UeHGMi6EIqS4cCbquhRyXJIQZkhD3cCVDMiunkOmJz5YfiUA69sOC09ljwo3\nlptcqalCOKxJT4blnKiBSGj8lQcsnNOj+CjRt2xaDFZXeR97Kv0yE/XLH/C6GS5sElgxjklg\n4CyPBQUK1AOeQ33USnFxIVjducLUxYbKwfJFGk921TdVRQa9mDPsE4QxmAGdCNOXbd2aAxOF\nGN/FS+Pl0nJDJVc3rCd/pBFedwVXyEifRxpW141UpEiTm6GoJ3xrI6Dx5K2YtB+WLZI4BPpM\nVCZ/pKFY6UaaKxi4Cd81oK6/wqduUCC0oiENkdo0jp8Tt2KKkowKY/V37y84vYzlZPWCJaz3\nTt4ag5lq3A6g4aGJuDMz3DONDWg4zlHMdNd4HaZsxTQtBWvCo3u/J7Svhi9l3hVY/7zuCp3X\nkSQ5lxJzRP2UPUlKzJkT51BAo/WusVC4PVhTHt17E82zP/2z46fYu4fExopwpWMZ+BEFird4\nJfEXE/bQ6SLHHQFj+jUqed8erCm/zmnQzbHew4+xqBcsOzvxuitUvBKOwxKGL4KWr7gl9fIQ\nPPUX71ppsuToHt48x5oycBDxC9YMFh5P/QO/RGbCdIHKCeFRoW+cSMDCbkz7PLJW8B414Nt8\nVDjl0b3H8FjYX+AsqpujXo0X0yEOJe+e0SANiwY/ByzaHV0n/O46tMGtmLkXFi4y6dG9zasv\neu04x0JCUZAkQP1UcGMSgphAMZFUFtS7z48BaCb7F8Ym2JDutEmaCtoCsBzXPKYzan0a9HQ7\nACgws6AAABRfSURBVIskVzpHlxJKpXg3mhc0cnODn+PKCFukK+mAtSy2TV98TbDk362rYzWX\n+27rWFo0uUI5D2RXME8aX2PAoq7ID5aOhXhZOBJO/r4ktk13eEtDYR7VDhapUeGhH8ygXJnj\nLjQqkkAlKV66woBZ67kKOKyFEtb7hEVmWMmousHCyZUpJnSf8Cln1ZiCZcdAZ7jnhcsZFS4c\n0TlaHyyz/s3o41XlIU7pSLsiig625gqwkAgskrkP1K8gDiL0JJ2YFqstwfpEbMohTunYhXaV\n+LRH/B/2V5ANm1TJhDtaPLC9likqAFASKhbmNa1WzbFe5C9o/MZrhzil43IlVeJk3JWERJqi\nJDAbEBalOcXj+i8HLA1yYq06KpTYA53Gy1kHKJD63JXCR3OFx3y64mR8jOLIHhdaHLnCY0sp\nx/ffDJ82dZFUOVbMcsMLWnuqQjlYQQSBwSA4KOdgGecU4geWceIefh07KgtrWnFaCNZlwq++\n9u6xXKxMRUlnV9gj4Ra6HY56XmdFoqSEqQYsD6+OpidMM7Six9r3KZ2Qu4J5wIVyThQHnSS5\nKbsbCsk75gldQhMUGeJliw4LwTpNOAm951M6HqyENc/yKbTihL45o8JRQcKOvVdQCKyMQXEh\nWN8JJ6F3fErHwcrvrqSEsrpJ3k1gQgm8rkN4IDJeCRccJHQz7oIwWNa6JtTiUIj+btKpLrA8\n3srmqpuqgxW0QBxoPtDuNCThd0OXCYcSuIrbc4DyhIWmisFaKj9Wan+QSxwQLIYcaYIYmjUc\n93QaBTFUTgQLLxi70FQtBCuTKgLLEwV1Km6lXvTMscFD6jQe0iopfZRJB6yZHgsPFyYsNE0M\n1iK5WKGjZbL20FE0g0HwWiTeSZOOoxmSkAapmd15jMrMsQT81Rw0FPqwMjsC1dpDcUeh44KF\nvJpTYnC/KcPSQBatMkeFU8Fyk4WUa7WBXKxwRqz9lT3dBsvk7SRzN0EOwp9HsKikQ4ZoJfYI\nqONVFun02BdYPqzc7Gp4bA+TgD8AiXyGuQo5xF+JWjPHejexj3IqdGdh+d1V9yZQ1o7BcuMO\nTIKs3A1/JlSSAZ3rsopSArD+zqK5RRXg38MncsImilPQXfUfIfkiAdADAUyyQx6qfGKwwGOZ\nAn3aTUukJWC9f0Q9frS0aqLIeoi4WyyXua+MPFjhVbaKDNbcMSHAoDZBXBU0kSR7K0oLwPrr\nNvp2bt7ye471RdnWakUNuyvn0mQ5xau46bkVL5Uvw0lYkXtrAVgdTDch2hv6fYcvg1lhrVbT\niLsaPsEzKvBWkngsBBKEQ0jqy9PCcoP5YzxKHcuLVdBdWc1wCPNLB06LIFyZkPCdzJqnbLk/\ngzVJPmgG3RW0gkHewN4y9SwaFKXtteza6kyVWyA9Flhj7mqAK6muahhJ5jFYOIGHuTT8LXQ4\nk0cWU7uetcjRwPJCg4/sQBhERNGDaZNhPBGpZUm6DHxZtqOE9Z5Si8Ci45dt1yq7/FhJ+Jsa\nwCoElhuKcK2BRE9pg5UgijFY2yvgrvo30yLIVQgs9AqdQoal0ypJ2qUJg26HSbUArIwqDazB\nKCjIOZyQTI6lx3X9K87EoTv7v+ljCKcZoJWZY2VUYWB5kRH6QLf/2zupDWFFRoWCIKnHgFLn\nVqTWQIIBXdruX79MUJmjwowqCqxwFDT/x9yVWgSBBJ0oNvtP4Nb8SQb2dxL1Zd6nbVqKeBrq\neJVFCjQRKz8yJusRciy78iwGWZJ+wzNMjcEdb/uWV0NIPUiPWo8VxGANy48VKRb0F4pG9eZN\nv4UDls6+nJzKXl69ChOYy9lzDNaQ/J6IxiZ8oajd0I4zLhi6FopmoP92ak+XDyxUhhissAJR\nUNht/vlX2JcZW5kQfMXOh/z00DgxkoqRbrSrEqXsuFYMVlDj7gq48oLlaW/BJoT9AicDTSjE\nXZHlKVhLzxsmFoMVUAArl6vQwMoKe9CDk2DBd/Ift4EyA7aG3Jw1Z3sxWF4FxnkiqhVt698W\nb5wEr0UHCKGrkP0LlCEGy6codzXM1RhYnnk4p8JuyeT3Ljt2HascMViuAryMYWUfXH+OZQ3v\nrBM2niV06YEUS8sDyRaDZSuE1QhXbjjyTdFRTXiaDIKlHZeV7pcrBstSnLvy+CtPM9eHCfWK\nwfJWD+gyxl3pHnyrVJYYLKKZ7moko1KA6VBmyLJ8lZ8Xk2OZKEqGlaWKwUIKXq/utgv8Ctq3\n5s4ZPfqbLljIH+HsUSFAWVTdyhaDBYp0V4O/qPCBpV4FnNRTIz9nodAVMXDJDYCl/xUqBktr\ngbuCZp7BP6IHn7FR08ylDH7bqECKk3eAqliyGKxe4V9tuQ3Dfs03ZsNg2RfqoVq6By2gyro+\nC4VBBqs0E0ThXwN6WoY6oRfb+cAyzgfmQe3BA5Zuo+Mm1MDgxFBobbYWgyUTuCuykPWO+ZB0\n08gcGg7RdX4k+UdDS7u3ssRghZ3QFHflLmXHO1/BVL3qzAsV1vFymip6wUMgfpaiw4MVxmoe\nV77xYeiMHgaLDPJ0QUGn+hAM0ZpxHatIE0qT3NWUC9sjNgElZXSQp/MnYQ0lTRAsmSilY4OV\n3F2ZZce9icnGhbAHeRgsfDGW6ZnBmql1dlyYFdd+PFcyMkjheGYN8iDVwnUsmMBgWSrrYePT\n3NUErqJkpWLCM1G4yVnpY0HQimCV9bDxbO4qUjZYJH6CG1Nx0ORXhY8FQSuCVdLDxgewWocr\nX1XCqmNJ6q2Er1m5WhGsgh7du7W70paMPxpaFU/9Qk8pGLIVwbL+InOYiNP27qq3RUd7boPA\nO3Qgx/bjhjqexxq8O9qU1oul6gx+wxFgeaeWonVzrAIeNl6Gu+rsWf+ltEaBZmxomtuL+6YX\nojXLDQU8bLwcd+UByz6pGD7JiNeXwdr+YeOD17xMap5CLlhqMnit0ElGaM1gFWBiUhTMjZXU\nRFlcgfcau6Sdc6zIbslFkuk1LQquwBU+W4jXQn8RI7uaR4Van6to7lI+TqIZeaRT+r01MQqm\n4Cri70OHO7gkRnqT+vkWNtOKYH2b1hc97luc0om4eQdtngArGeNPsNdCOZaQNDzWp1XLDe3T\nwhpx/crvbdVywxBWucJgZAaE8yxyDXI1lzEEtGqBtFtadIWGFQukMfcaim0fLbWtYy7LCnv2\nKLBertY/pWN2XQ4THg1ilS9r1yn5CFnBfKrszDxCG3is9nXkwZnpdmjMHaysBZKMBk2SNBOs\nsjPzCG2QY7UPJl8px9rIXUH3YyUDu5a1H+15VBh1w734JaYp8vp0q5a1H+24jjWIVW6uIHMf\n25bBa7LqVTmV98QmtnJXwIgg3+y5w13UD9pOwRrGKh9XeDCHruTzzB3owl6sSu0TrJEnB/oX\nSVRlQBbwZVXu3GAXkNBXrD2CNYxV5iIDehfeqSNbR87sVKz9gTX2oNPpy8SLwTLaHVgz3NXY\nw1HjxWAZ7QysYUTyuiuw4OdreY5VU0a/K7BGomB+rvw/aA7MDXYRGBXWdfpwT2DNcVfpwqA2\n4/ygOTw33IWvWZS/K0b7AWsEqxXcVTqhKiupt+L3wrUXsEYACbqrErmCmEfqrXruRms1UTsB\na5a7Sh4GEwlAEtZ3/F64dgFWae4qevQ2kEvhgaEpQAj6uOiStQOwxrBam6vo0VugYRAs1byO\n0WH9YM1zVxnDYPToLdBwzGPVMTqsHazS3NWEXCjYMC7HKpysusEaw2N9d5UErKhRIYM1Q5Em\nZrqrvFxpq/PBiqtjMVgzFGVivrvKeS81zcL8HGusOedY8xVhYhSrLdwVFAMWjArHmvOocL5G\nTYzSEewhKVfOqUD/5OjlI5tzHWuuxkwscFfZLmXopljvB1aNYBXirjzJDoNlVB9Y41itypWP\nLOaqQrCWuKu0absXrCoS6zVUGVjFuCsZins1JNZrqCqw4h/IPGfZyeK4N6CawBrHarUwCObY\nQflVD1iFuaveImMVUi1gRWC1Plf+1XBP8x1RlYA1TkYpWLkvx1QlYI0uUAZXgUupjqh9gLWk\n7pVSUIE4fA1+D2CV4q4YLKQdgFWKu5IMFlL1YA26q9V/N4hyrNjL/Xaq2sEqyF21govxdnG/\nxwWqG6yhA7fRz5zBUx0Zq5XB+rt3T1gVl1uaJ6wOYrV2eoUxOnyGtSpY3xN60mWCBwiU5K6s\nciiDtSZYN9E8392nJE+xX3QxYGpZ5VAGa02wGvE2n99LH9JUkrvygHT0uvuqYFm3ultkoih3\n5QPr4GcKK/VYg0dskyqD9S6PPiZcO8d6fbpPC3Os0rCSHPpcrVluOKNR4ek710Rx7qoVhz5b\n69axbl0dq7nc59exBq1veO/Ho4c+W5VV3ot0VyyPygFLYM0xzVyVpDXB+t7aoeD9JMT5OcfE\ncLApnqtjBcsVwfo0vz3bPxh61imd2rGSh0rvVwTrKi7f38v10z0demq5oXJ3dbiCxKqV9696\n+UXFqQXSYaPFY3W8s4drn9JpBPoSbaJ2d8VgZVqk07U9pXPvz+t8h5MsMfjVUg1cMVh5Fun0\nFs3tLS/Nj6zXSbyiTYxkvDVgJTnHyrNIr1cDhap7tIkduKtWPCrMsYjW89pdRXq5f2JNjLmr\nWriSXMfKschsEyO2ct6w/UgUZFDRYG3mro4WtzKoZLDG3FVGfxVjnzWkcsEacxglPGeJFVSx\nYG3nrhisFCoUrFF3lXU0yGAtV6FgjczPXWTgHGuxagQrf/GKR4WLVSFYq9REuY61UPWBVVGt\n/ciqDayazuEcWpWBxVjVorrAYq6qUU1gcRisSBWBxVjVpHrAYq6qUi1gcRisTJWAxVjVpkLB\nYlWvGUc9PUgF2V3Fyo42JaEVBqsOI9VZYbDqMFKdFQarDiPVWWGw6jBSnRUGqw4j1VlhsOow\nUp0VBqsOI9VZYbDqMFKdFQarDiPVWWGw6jBSnRUGqw4j1Vnh30WxsojBYmURg8XKIgaLlUUM\nFiuLGCxWFjFYrCxisFhZxGCxsojBYmURg8XKIgaLlUUMFiuLGCxWFjFYrCxisFhZtCZYD2UM\n32ji1ojm9k1o42T6Q10ntgJGMm7K9yrE9S3trvNZSbstK4L1Viv9Rltw7j6dktm4df01X9p1\nYitgJOem9E+wfVtd57OSeFvWA+vdGLAuetqfaN7tjL9UNsT123rGK+k6sRVkJOOm3Nr+b13/\n+TYFW0m8LauB9RBnBdYDniF96x5T/hx5qHS8Lr2F1hDqOrEVZCTjpjTiq4xk3BRsJfG2rAaW\nuEkD1kNPvIj2edLobyWRLUG6zmOlByvzpohGrrEprZXE27IaWG+pwbqI1/WXGnbW4Y8/ob7i\nTLrOYqUzkntTbt3Bzr0pvZXE27LmqNCA1eksc4H1aD157qPRGcm7KU8h7OOc0UribdkALCGe\nvz/47s8kyyH/NBeZHSxtJOemPC5Nl+bk3RSwknRbNgCr17cdzeY45N/mjPrMdDSUEfUl16ZI\nebWOc0YrvZJty2Zgdd+aDPvp3JdfUNcZrJxJjSfXprTHucm9KdqKUior24LVDz8+CQc5n9P5\n031AXSe3Yowo5dkUX9cZrfgNzu9y8UpNsNUb62sn3Wrfu4LJq88eU+gldIhCXae2AkYyboru\n+pR1U5CVxNuyAVi3doW/XREudSH5Yw55xnI1MpJxU7qa+PfSZj+5K++9lcTbsgFY3/78VPfn\ncDJD3CS6oicpoK7TWkFGMm6KOotnr38+K4m3ZYsc63trxOlhPjbJAiE8XYV2ndaKbSTPpnRX\nGLhdZ7aSblv4eixWFjFYrCxisFhZxGCxsojBYmURg8XKIgaLlUUMFiuLGCxWFjFYrCxisFhZ\nxGCxsojBYmURg8XKIgaLlUUMFiuLGCxWFjFYrCxisFhZxGCxsojBYmURg8XKIgaLlUUMFiuL\nGCxWFjFYrCxisFhZxGCxsojBYmURg8XKIgaLlUUMFiuLGCxWFh0brO+jfR7D5THQ5OWdEHOn\naqdN4ptol60jbaujV6Nu+th8Qk1O9g7qJzBYYzrSttp6CXFtbwz8d0E30Le0gA4G66hqhI5z\n6JkflhismTrStlp6wpMXPrfulubts557wl5nIc4vdYfk9vsvF+tuI6wm9Iz82ve3GRbi8/N6\n3XMjTUvC0a3Rz2vU87/qybjduza3Ix0YrIt9g/yLubv5o0+9Hpqju9D3QMdgndFz2Lps7Y5b\nYrC6lhfakzLfPsfUmNuRDgyWHZle4vyV33MbH5v28dtPeBBW98i1Z/dZmOT9qZ7f8Gy//hZ9\n9O2tlq10S4Hnv9qHQrRR+IPM7UcMlnkawKV7mMy3DZDCZF+YPgusi3riTPc41z/S1gKrd04v\na/6pfx7zCZvbjxgsAxZ64sTtF7jeb9RIfl73swWW8xjB/tVqSSyR+d3Tvf/6R4hrc/vRgcHC\nOZYFlrw3qrqlmDibOWNg2S2lpGDB/O4pgffuCW7G3H50YLCefZLTCYGi9LqdIMe6/kZ/r08M\nWE5LKUlLNF+2T9o6nYi5/ejAYKE61hdyJiTLF9lg6Rzr4tBmg9W3/LPmy7c4v7vHMYO5/WhX\nGzNRv2z60lXefzlOo8Zuv8zn0ubVTzMqbANUm5y/debUT7BGhVJqsFBLMIRGhWZ+a6bpegNz\n+9GRwZJ/+lxhHxTP5rzhs5/61z25r+myazrBqWNJ9YpaIg/UVciu1vyO7A6mJ5q2Fx0arN8h\nvfzYOt9V2vz4UXPtPnel8M6bnbrziNf2axf1+gmq8t6Yyrt5hZY4tN1N5R3mdxG4j4TG3H50\ncLC21auPq7sUg7Whzvs6i0PEYG0mkfLZzsWJwdpMDVxdsUMxWKwsYrBYWcRgsbKIwWJlEYPF\nyiIGi5VFDBYrixgsVhYxWKwsYrBYWcRgsbKIwWJlEYPFyiIGi5VFDBYrixgsVhYxWKwsYrBY\nWcRgsbKIwWJlEYPFyiIGi5VFDBYrixgsVhYxWKwsYrBYWcRgsbLoP9NzMcLDajUdAAAAAElF\nTkSuQmCC",
      "text/plain": [
       "Plot with title \"95% Confidence intervals\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "options(repr.plot.height=5, repr.plot.width=5)\n",
    "\n",
    "# Plot the fitted regression line (the fitted values)\n",
    "plot(data$Gestational.Days, data$Birth.Weight, xlab=\"Gestational days\", ylab=\"Birthweight (oz)\", main=\"95% Confidence intervals\")\n",
    "abline(model1)\n",
    "\n",
    "# Add the confidence intervals for the fitted regression line\n",
    "conf_interval<-predict(model1, newdata=data, interval=\"confidence\", level=0.95)\n",
    "lines(data$Gestational.Days, conf_interval[,2], col=\"blue\")\n",
    "lines(data$Gestational.Days, conf_interval[,3], col=\"blue\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how the confidence interval around the fitted line is narrowest in the centre of the x-axis, where most of our data are concentrated, and widest at the extremes. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12.8.2 Prediction intervals\n",
    "\n",
    "The 95% confidence interval around the fitted line describes our certainty about where the fitted line is (i.e. where the expected value is). \n",
    "\n",
    "Sometimes, we are not interested in the average outcome at a particular point, but the likely spread of values arount the average. In this case, we are interested in obtaining what is called a **prediction interval**, or **reference range**. \n",
    "\n",
    "A 95% prediction interval, or 95% reference range, is an interval within which 95% of future observations are expected to lie. \n",
    "\n",
    "The predicted value for an individual with $X=x$ is the fitted value, as above. However, there are now two sourcces of uncertainty to take into account. (1) There is uncertainty about the fitted value (the expected value), as above. (2) There is random error around that point ($\\sigma^2$). Thus, the variance in the individual prediction is given by:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "V(\\hat{y_x}) + \\sigma^2 &= \\sigma^2 \\left(\\frac{1}{n}+\\frac{(x-\\bar{x})^2}{SS_{xx}}\\right)+ \\sigma^2 \\\\\n",
    " & = \\sigma^2\\left(1+\\frac{1}{n}+\\frac{(x-\\bar{x})^2}{S_{xx}}\\right)\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "\n",
    "A 95\\% prediction interval is then given by: \n",
    "\n",
    "$$\n",
    "\\hat{y_x} \\pm t_{n-2, 0.975} \\hat{\\sigma} \\sqrt{1+ \\frac{1}{n}+ \\frac{(x-\\bar{x})^2}{S_{xx}}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "The R code below calculates a 95\\% prediction interval for the birthweight of babies who are born at 280 days' gestation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>fit</th><th scope=col>lwr</th><th scope=col>upr</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>119.8818</td><td>87.01496</td><td>152.7486</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lll}\n",
       " fit & lwr & upr\\\\\n",
       "\\hline\n",
       "\t 119.8818 & 87.01496 & 152.7486\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| fit | lwr | upr |\n",
       "|---|---|---|\n",
       "| 119.8818 | 87.01496 | 152.7486 |\n",
       "\n"
      ],
      "text/plain": [
       "  fit      lwr      upr     \n",
       "1 119.8818 87.01496 152.7486"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Prediction interval\n",
    "predict(model1, newdata=new.data, interval=\"prediction\", level=0.95)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 95% prediction interval for babies born at 280 days' gestation is (87.0, 152.7). This means that we would expect 95% of babies born after 280 gestational days to weigh between 87 and 152.7 ounces. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12.8.3 Comparing intervals\n",
    "\n",
    "The code below produces two scatterplots of gestational days against birthweight with the linear regression line of best fit (obtained from Model 1) superimposed. The blue lines on the left-hand side plot represent the 95% confidence intervals for the fitted values across the entire range of gestational days. The blue lines on the right-hand side plot represent the 95% prediction intervals. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8AAAAJYCAMAAACaSn8zAAAAM1BMVEUAAAAAAP9NTU1oaGh8\nfHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////UNI3wAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAgAElEQVR4nO1diaKjKgyldp/etvz/105VloCoLAFBc96b3roAgePBENAy\nTiAQmgXb2gACgRAPEjCB0DBIwARCwyABEwgNgwRMIDQMEjCB0DBIwARCwyABEwgNgwRMIDQM\nEjCB0DBIwARCwyABEwgNgwRMIDQMEjCB0DBIwARCwyABEwgNgwRMIDQMEjCB0DBIwARCwyAB\nEwgNgwRMIDQMEjCB0DBIwARCwyABEwgNgwRMIDQMEjCB0DBIwARCwyABEwgNgwRMIDQMEjCB\n0DBIwARCwyABEwgNgwRMIDQMEjCB0DBIwARCwyABEwgNgwRMIDQMEjCB0DBIwARCwyABEwgN\ngwRMIDSMDQX8vjF2fo7fv0yg33iw7iH2dnaiv9uZse76z6uER8fY7VdHZtTS2ozCXB431Nx2\ngXCa5Und9eWRv8jN3YY3Pn9sMbuZnEKRn9ntrpzXyNJl2PgDzA4HBuru7Gkluipy/9ZLePQn\nlhTwXxeX854FHEEz07ivF7AgYMlHuoDrZXazK+ctObr2Ww/A7JX9/Y17O/Y1E10At+/VIs7i\npGICjs15xwKOoRmQzNb76QUBy33pAq6X2c2unJ9jdf3ytxDZBVDV13mo99PugH/33+75I/vz\nlFfEImZaL1+j1kvzZoihWTbHaxwBrWCh8fDatV5mN7tyfuT0/e6vh+4HQrCiitkz+xhJfud2\nYs/7/Bi/vW69m/ySKRl//S6S21ttwr+fW9ePxuTm996x7v5xpfwd+11ylxefnKhsZJNk+uZi\nZ/05s/ufvBhv40X8rx8NnMdzZKpH72F4ju8bQQTN+qz3+E20ILeYMOnUrErmNB8qw6WrxSy8\nGWY3E7BW13Vk6iabUjL7sm+ydzYZE0uf+irzvCv/esKfGIBdxOanA06akVIfu09ONK2fKdDO\nur+k+ktZJO2A4cM5ZioxXtwHImgGMlf8DS1otqtFp8Xq3SXghavFLrwZZrcUcN81f4cqPmWl\n+85RDo4u9gDo1zRWX61iWiMnakuErkz+On0YbnaTlODUl32itn6hQEfW//rup8/tNVxcvwpf\nvkOXdFW53fqz+Pcy7aYaRgTNxh2447oFzXa16LRJfk0vgKWrxS68GWa3HAP/KvodO1DdtB8V\nnnyzs5VEESvRn/kbE38f4proW/Y1ZG14VePff+PBfmA1jryG8se2tVL+jnXvwbazfaJpirvA\nSdYX4UXeuPSzzqIvghfYWPv+Wrcr3jAiaJ6MgWULGu1q0ckdzFkXwMrVYhbeDLPbR6H7Kt4u\nfdO+L6PT+uj6AdON2bOAEwHfZJd2V0z3ab5OAV/FpMVLbX7Hg9dJSnHqdxhomyeaprgLdGY9\nUPu1SYQ095fizWfmsyVE0AzufqN3K1vQaFeLTu5gzroAVq4WXTifHqyX2e3Cn2KCsAPN9wZt\nMMzuPzox1z9gImAm2pN/YHMZgQj9V6VWmwLd/Kl8cqJpirvASdajkY++3/43RnN6m//dLwyW\nJ+ZYdqbhcJqhgF9i+2vu7yZ0upibsrp0tXBHdg0wu52AhyU61/fZoYse/ey+nusfcLYnf/Xp\nKQKeppwRsLlzWsKUZiOzTx/EkOP4f+fpKXd5cVgj/bYRTLNqv8v9a5xutGuUgI39SQKuhtkN\nBSwM6L2R23XSlv3svp7rHzCJQjPQp07vo8sC7ryonpxonu8u0J31EISTN59fd83Ot+fbLO/7\nbwxh7ikMLRBA88TTkttGu0YJeOlq4Y7sGmB2MwFfr0NzvgZRXsYe+E/X8DlOBIz/BP50F/Yn\nh6f2qGbYdJIgB03/zE3Op6derDHw1PVZpNmd9VCytPisvUNuXkIvK6jSOCJonhWw0a4WnS7m\nphfA0tViFtYMs9sJeGjF1zjRPwYP352O9A6z+zazaiVWP6i4OuOKw3lOEp5jXPGfCFv+G9dT\n/xsvpumpKpZpnigwT/N3Lmu5lP8LEhj99FkdmzzC0S5iaJ4TsNGuFp0u5iQf4tjK1WIW1gyz\nmwlYBDfEyF/Or53V0THMZzL7BXN/w71Yr40epvIWBTw3Dwwm3KenPu0TBeZLuM9lzcdJibEX\nGAOxcg5EXWCXzxDx8FjC3wpiaJ4TsNmuK/PAkrm7PrZ8tZiFNcPsdt6amBQcxz5/QJU9zkML\n2YOjcT5xvATGM+W2+diYmwRRxlVsykvr7joVrsQyThRwlyBpdGYt9v+D9e2Ds3/qFBnq2NUQ\nOILmWQEb7WrR6WJO8iEzWLxazMKaYXbD4dbz51moJz77ha3dXT6VImb37fBkv+t3HjvrgPyw\nDVe3gr82Q59+VQFcC60NsE/tl7wq4+CJ3CjJTnYVF4cra/FdVLGPzXa392e4tOUpwyjpsqN1\nWD3CaZ4VsMmESaeTuas52l28WszCmmF2QwEv4Cr4HOf6CTsF0ZyOOgVMIBC8QAImEBoGCZhA\naBgkYAKhYZCACYSGQQImEBoGCZhAaBgkYAKhYZCACYSGQQImEBoGCZhAaBgkYAKhYZCACYSG\nQQImEBoGCZhAaBgkYAKhYZCACYSGQQImEBoGCZhAaBgkYAKhYZCACYSGQQImEBoGCZhAaBgk\nYAKhYZCACYSGES/gv8f4szfX+9/6yYTmQPw2gVgBf9XvkO/t17gIPYjfRhAr4Dvr/r2Hb59X\nt6ffwyQMIH4bQayAO/ZW3997+kVqwgDitxHECtj6SUYMUwgVgfhtBHQHJrhA/DaChDHwa/yd\ndRoj7RHEbyOIdo4uIEp5/q6fT2gLxG8bSJgHvg/zhN31QfOEewTx2wQoPEEgNAwSMIHQMKIF\n/Lmx7sH588w6CnHsEMRvG4heStn1A6Tng5ba7RPEbyOIn0b69cv3jt2+/HtfnmZgBDREshUO\n4jcHTmsnhBMVv5Bj5G6YYFie6KdhNhrKNSXxmwGntRMKCpgx/elYapfYrRDcKNeUxC86Tqv6\nLSngDhD8pR66EMrfgftP4hcDHvotKWA5Rrp/xXf8IggTlB8DE79I8NFvSQEHRCmJYDSUa0ri\nFxde+i0p4IB5QiIYDQWbkvjFhJ9+iwq4qiKOgiqbskqjqsLJU78k4L2jyqas0qia4CtfEvDu\nUWVTVmlURfDXLwl4MxSaDa2yKas0ChkJ/AbolwS8EQZ2S0i4yqas0ihUpPAbot+yK7G8F+Mc\ngGDwWaCgEiB+ARL4DdJvSQE/iWAFZv3NX1J+EL8aCfyG6beoC/3ufB8yI4LRSyoA4lchnt9A\n/ZYdA79931VIBKOXVALEr0Q0v6H6LRzEeoJXB2cqolIMPiUToQ22yzEwJ36T+Q3WL0Whi0DH\nJBXHFIXeEZD49V5+BYsOTkECDgd40E51zTQPvB/g8BujXxJwCUhmh/8lwyTg3QCH3yj9koBL\nwEUwudD7AQq/cfolAZeAs4fWn/nLrg1VGhUPDH7j5EsCLgPXGMn4m7no2lClUQlI5zdWvyTg\nIphGKeWB/EVnLyECVRqVgGR+o/VLAi4Ea56QBLwzpPEbr18S8EagMfC+EcRvgn5JwBuBotD7\nRgi/KfolAW8Ga54w17RwlU1ZpVHI8OY3Sb8k4DqQ74ZcZVNWaVROzPMbOf2rcy6SpMIi6kK+\nIXGVTVmlUTkxy2+qfknA24LBycMs9a6yKas0KgfW+E3WLwl4SyjHigS8S6zym65fEvCWUI4V\nCXiXWOM3Xb4k4C0BaKUx8A6xxi+GfknAGwISTFHo/WGFXxT9koA3hOFY0Tzw7rDML45+ScDl\nwayl7nkrW2VTVmkUFvz4RdIvCbg0oDNVYj1llU1ZpVE48OMXIfwsyyuSpMIitoLZK+d/r06V\nTVmlUTjw4hdPvyTgwnBMKOQVcZVNWaVRKPDiF1G/JODCmBCc242usimrNAoFPvxi6pcEXBhT\ngs3NbAVWhSqNQoEHv6j6JQEXg7Es1uI3Y42rbMoqjUqEL7+4+iUBF4JeFmu6VCTgfcCbX2T9\nkoALAXTMRlCDBLwP+PKLLN/CAv57XIefjr3e/3IVUSnAtAJzHdnHGJj4XeMXXb8lBfw9g59/\nXv4l2b0S7AhJ7icKTfyu8ouv35ICvrPu3/jrk59Xt/xLsrsm2K7dXuaBid81fjPot6SAO/Dj\nsW/W5SiiWhjv+zYVm3kxVrmmJH6X+cUOX4nsiyQZ07G5DbQiKgVj4+BIsszAq/sze9AFm5L4\nXeQ3j37pDlwCMjgpHSxmRyd3IWDid4nfTPotPAZ+fYZvBxsjGSFKSa6p3D0EsYjfBX5z6bfo\nNNIFRCnP3yxFVAkd4dCO1g4FTPzO85tNv4Xnge/DPGF3fRxqntCcY2Bw5/z8IXLhRUD8zvD7\n029F/BZo/X0RLMdIYMMeI+ULZVXZlFUaFY81fsf7bzX8koBDYejT+gFZOP2QpewcmaaiSqPi\nscKv9J9r4Tfaju+9D00+zoxd/mUqoh5MJgN1nyzHifKAOpij2gWbkvgdvtr8av+5En5jzfh0\nv3p8u0MstXP5xJpEEetgTA+bxm85LMmQpxvEr/w0+O31Wxe/sWbc2PX7+7h9flzfdj7N4Opx\nNenjL7ireKUitxKCI0H88im/w+23Mn5jzWDsKz5+3ta+J/pnZofUVL8MXYlOGrCezZT8IH6n\n/Bq331r4jRcw75frgA3rMEBkEbVghmC1VxE8fEonuhaCY0sifuVeyW8fvqqP33gX+s35Y1xv\n910eJB2FYCngbPSWdaGJX7FX8DuEn+vjN9aSN+vub37tfgy/zuyVo4gqIBbUcUc1VNRK/tGx\nkFxVLteUxK/F7+g/18dvtCmvTvtQjzxFbA+wHsfudNUhEKUEz7BksidTvg4QvxzyO8avKuQ3\nwZR/t+GtDdfHJ1sRW0P2zQ6fCXTbxhxS1kFh0aYkfjW/cvlVffwWaP12CV4YHSnqufSyBNfo\nRsCF81U2ZZVGecGb3yF+NezIW1kSMC68CJZByUzOlfHgS5VNWaVRXvDldxj/5nWeZw3JkKTC\nIjLBj2BuxTtQYT64VmVTVmmUFzz5Vcs3sk+ZkYCRMROfhIdgCBqIGQmWfqtsyiqN8oMPvyJ8\nlYdfq0wSMDIW3CbQJ6uZfXRH2r7/VtmUVRrlBw9+xdMLefi1iiy7Est7MU7DBC/2t6ZfJQPR\n4yGcwif+c7mmJH7HQ+rpoxz8WuUVdaGfByF4Cj1jpKeALX6R6jwd/5ZrSuKXifBzPn6tcuOy\njbbk3S0/ZIZQRI0wpv4VseoDk2Az/qzGYIVA/Er/ORe/03LLBrHeyw+ZYRRRIYALpV0r2Tuj\nEuyaPyrZlMSveHwhE7+OYgtHoZ/g1cGZiqgFyokEsWdzAaV6ahSclgTn/G/Rpjw8vzL+nIVf\ns3zrb0TSjGie4MlTCvKRFDBGVJ00UpTSvX6jyqas0qgQzPCr9JuDX9OAyZeItPnQPsH60xYw\nZ9zsoRfjmv6YWX9VZVNWaVQI3PyO8edM/DrKN79FJM6G1gk2/BtjjKS5ld0zUl3n1k9W2ZRV\nGhUAN79y+UYWfp3lcxJwHpgEu6LQRigaocTZ9c9VNmWVRgXAye8YvsrEr7N4+3tE8kzYFcFA\nvQzcg7Vnlc9/xskcH1UaFQAHv2L6KBO/ztJjMycBr8PNHAhwyG21P6m4heePqmzKKo0KwYRf\n8facTPw6ynZvRWSQBe0T7Io8mkFo7V8lhymXnh+ssimrNCoENmWmfrH5NYte3IzIIQeaJ9gZ\neVQulopUcs7VICm+zsbz+5O7fnS2GVGlUWEw+DXfvoHNr1msvSMiDxRLti6iPMzIxnSOIbbS\ny+/fqLIpqzQqHkK/mfiFQOGXBBwHF8EsneCV9+dU2ZRVGhUN8eNlmfiFwOGXBOwNI3qh1lGC\nxbLwjLhKr73/qsqmrNKoCDARfx43svBrFOe5LyIbZOyDYGZ0wDDGYax3nwlZ+2H1/XVVNmWV\nRgVj5Pd00pv4/Brlee+MyAcXeyFYf3LOLWrlZ1KUcv39k1U2ZZVGBWOgbHx5+7iJz69RXMDe\niIxQsQuC1dhHEwomGbj6wp0hay94vD+2yqas0qhQCP3m5NcoLmh3RE6Y2A/BXHXAcoGO7Kg5\nS+6Vfd7/XGVTVmlUKOTwNx+/RmmB+yOyQsReCNZ/FauKXT06ioXX+9urbMoqjQqF/vEFnodf\no7DgAxF54WEXBOsxsBnbUDzztHdGAv1aXT39MkMRKP3m4Rdg4VZOAs4HHctgM0ipKNTv7JFK\nm7JKo0Ihf7w7E78AS7mQgHNCBjJ0t2xHOqKzntXvydBvnU1ZpVGB0G+/ysIvwGImJODckBSr\nQCUX3pXytGKilPP6tYuPMDk7qjQqDHL5RiZ+AZaTkoDzAoyOgLcl5wdV+CN0nnBOvydbv3U2\nZZVGBQHKNwO/ACvpSMB5Id+FpaaDOeihjZUeIXX212+dTVmlUSHQy6+y8AuwlowEjAjV5eo9\nXHMJ/SnZNcO+mZlJjTxNBOi3zqas0qh1KH6BfnPwCw+v2hRUg9gkFRaBD6BeRYoiEgY2jFET\nINjharm8r1OIfutsyiqNWoPi9wSWT3J8fo3j61ZFVCQ8SYVF4INNPo0emqk3y3I5XLIJNpI6\nchswp1+nfCttyiqNWoPkYgg/g33I/E4Pe1gVBBKwBTn2Gb5PFAxecab/6t5YLbjTh3W28DoR\nCNRvnU1ZpVGzMPmVT++LY9j8GgX7GBdRn/AkFRaBBknVjICBd6U+9A6Dd/gBv8L2ADplPvqt\nsymrNGoGFr/q6X15FJdfu+B18yJqFJ6kwiLQYAvXpkTzBZfacQYZlB6XdMk4GHEZuUH9QiNm\n9VtnU1Zp1AxMXvXT++o4Jr+Tcj3Ni6hRDP4e16GK1/tfriIKYer9SIcIMMJAD8yMv8DL0tkY\ne7XfBU7i4qcrYbnj7gVLE2oZiv3za+gXn19QsG8TFRTw9wzqtvxLspUTDGOHTO2SHpQ4pjtm\nJ3SUUqQHa+ItgtVJY/TT5nxRvwWb8gD8yg40E7+gYO8WKijgO+v+jb8++Xl1y78kWzvBk0/o\nETFwK1Z7JG2KXX3PNtKzsbMGoUt9kpy9YAH6LdiU++dXdqCZ+AVF+jdQQQF34Mdj36zLUUQZ\nmN6VINI8BHtc278y4hwc9MK6W9dXCWgJ+PIWDnd7GZsfu+cXTB/l4Be66+HGRtQvOB2b20Ar\nogwsgoG/ZTnUq5AdsNgAHTp0yEdYs49gt5ex+bFzfse3b3C5D53fScFhxuZNMmCnPbThBmkd\nS8JMMq0d6kJXnbKk2syWc7h8wFS1p7H5sW9+jenfDPyCAkNuwIXHwK/P8G1PYyTXLtXputnV\nLMKx0JBOf1olgOsnSL9lx8D75dfWLza/YEeQfksKmF9AJc/fLEWUgeE1O3ZNCXbCJFj30Q6C\n5fIBPXoa93rYGl3LYOyYX/D2K56DX1Awq1fA/O8+zBN218d+5glduyYEy024fzwRXBNMDYom\nBKunx035eui3aFPult+T/PFucZQj8ztbsIelYadHJqmwiKxgMlIx7ahFNyyHQHKMNEkvPxlj\nUL/6HC/5VtqUVRo1j0lXicyvfSAEJOAcWIlSmv03nw57wGwD50n6rbMpqzRqFo7XnKDya97a\nA7GJgFfdhCYIXvR23MRO+J3NRx7jofq1c9uiKffFr/sxazx+F13qWavUtl9lEpNMjMheRHY4\ne1bjBNVPAz7Bs6OrGXCtX9sPmx3+TjMlAcdBteT8c5o4/Oq/Hq2Cwm9s60/8DPwiCmK1zRe6\naPWstwfBev0A1K+/VeWacp/8Lrk6GPzqv173XyOZb6qZUkPx1+2I4InXM63RLL3Kxxq6cXdK\n8eWk5x/99Wu0X7mm3CO/zvC/3pPMLyjNX7+p/Ea3/vfKLsNMv5Ndb/argNWUju7WDnS4u2yX\nU6T3nNhUv0uzR5sKeIf8gpc/8xz8gjuqV5NsLGDO/zH2j+9ijGQL2Ng5bjA5mQAI5vq7nGeY\npNR7TlytH/DS78YC3h2/+u11efgFYvZrkc0FzD8Xdv02TTCb9J3ARzIbdnSh3IDMmykls2zU\nr+7Ge5xWVm9sOQbusSt+pX5z8au/eLskKPwaSf7u/QK6y9orGDQerHu1S7B2f4ypvClNIoyh\nQ5VzkOFKmHDM86ReP+x3+zXtg9kl4MD89vrNya9uJf/mwI5C/9MvYTi/PJO/z6xhguHnWI2R\nHifBapZBTzfYPpfsfi2CufKf4TGfxZN24yY15ZH5HeWbj9+FJdGLFqbzq5L8/KXL892vWv/+\nPS4igOGBW7MEz41BmN3Ncg4HSPNTDuIq4DBWOTJ7Ghy48ZiA3+LnGZMjcGh+jd4zA79Cw5Df\nJJMjkrzYHT5y8rkz307au4jaMDMWUgEM41wfcN19K9do2DiN6+f13kj9JjTlofmFq99y8GuN\nxdJNjkhytZ8Y+96STHEUURtmBDz8nRKxQCr4Dnphlb2Y/4We1lr4as3kcByYX9DWefjV2xve\ngTOiVoJdRMiREmCYSYdruXe2vxsXkNIvlx51isWVoUqjeqjOctzKxS/Ub2JrkICDYLpCXI+O\nZEQSnLPI7oRqDjp5oV/V+Sfot86mrNKoHgx0lvn4BVMYw1aSxWlJ1AuAPTx5u1a4VpUCYEF8\njhQBZwvsn+2ljWilCmhrb1nsF1nF6zfVQTsgv+aPL2Tgd/wK+U2xNy0Jkwx7EPzcBcEKWsdq\npKNJ5haxk6sbLNxR3bxoRD3aRdBvsoAPx689eYfOrwaCftMFfBsZ9omlvbvl9/UnWVUcgGAd\npdQEuztmm3JIs7hGTmAENuyNDV+ZZkamPhy/k4evsflVUHuT7E0WML+wm68Z7+V3FSZZVRzQ\nk5Y9rEXwLMvygPbE5BbQr8glTb/JAj4Yv0b4Wf1F5Jcb+bNE+WII+Mfw3bcfeYJXByNbVR7A\nAVKjGZO8tUCHJBfkZeo3fvbIsjI69cH4nf76MjK/Vu7pQBDwwHByTzJbRLUADhADcUXQ3/pC\nRSnx9Ysg4APxazR2Dn515lgmYwiYd+x+DIItMFlr4WYJn0l+XaEYBipFXvj6xRDwUfidvnwS\nm1+dM5rRKAL+dOwIBC9ABRkVv4DAeZ7BJdBnkEG/KAI+Br+rb0pI5tfMDgeJAhboGU63ZbGI\nuqE9LOUurY6QAMU8k//MkZryCPwuvykBg1+dG6LZOALGRo0Er0DyBAgWgx8PivsMTpZ+Tyj6\nrbMpKzRq/UnrRH51nbfuC80k//oHvq//kMxxFlEabPF2s3BUj3UUu3yFZXllMJd+w41z25SC\ng/ALpo/y8KuzXa58AX6NJPIHrXyn8COKKIyJw+N/1IhSrnfO8LCnfpeNc9oUcK4DB+EXvHzS\ncRQmjeRX57pY9yL8wiRP1vUPib469gzPyK+I0lhu5hUSIHky4DH+c5ANroPTSbyAyXanw4qf\nTRGLg/BrL99A53dMy1fFWYRfmOQsZu7f7ByekV8RhbHs6MwdZaazxMDX1VAHV/dfeUWs6Deo\nfdKa8hj8updfGanS+BVrOcBzDf7GedYnADCJduy3HpljFx0iYD3RJwn1CVBqyPvveFkshK+K\nC/gI/MLGzsMvEytBBL/Lg+xp8V71CYD7DtyFZ+RXRGF4CJjZvaj2kDjooz3mCYF+x3yEfvns\n7/EsGbeYJA4H4Pc0WT6JzS+TGa4+P1hewAcZI8GjsiO2Ukg/KZzaqX5DfvvIrzqx2D+/Vltn\n4FfKV+l3VcEFx8BHiVKCo0ySqfeJZJJgzqWj5UMwjF/J4e9CEGXJOGeKgHMd2Du/ztWTqPwq\nBTMfAZfg10zy78p2Pk9obsr+FOyb9NABQySpWHE9qGiWr3HrlQk6e4p98zu+uignv/DTx0Uu\nwG+B1t+WYAN2n+jiwBgjLXTOjn0n8f6cMcW6fsPtR8wLDbUYNXl1Lzq/4KC+BW/M78EEDD7h\nN7PX5sr5CgtTCv2a07+oEd+KmlKjEqN0Z5mLX63sMc9wF3kNCQI+xHuDIZ9sbH5jvzwcNCBy\n6FeEr/iJ4U7YJDTl7vmFzk4efrV+Va4V8CuTHOLN/ZpgFYhY6EQj5Cv8ZzEYRnh6YaYC4dg5\nvyf108uc5+HXkLAoBR0pLnTsb+dksSoTAMFqe74T1fP1i2Ta+u3TqWB0pgpEYNf8qqm7ERn4\nNeRbpYDjfr0usIitIVtehjAshiHZTPffK4xC/Y7XTC79pjXlfvlVzk42foF61/qFFKQGsYJ/\nPza8iG0BpvAUeXK34W4xuaLdz9MC91+l3xwUJ+a4U37h04M8C79AwRzMJaHXhKLQq2AytsH0\nltzgxieTYYr1J71Pav0zBxsZKl5VU0psbZT18roM/CrxygiWKgUZJGBfGDMBnMP+VF4ATDK8\nyi2IX8kNnRWu2dgZYmBjo+Yftcbi16FjnRUqSMC+MKKUXEYzxu9csuuJWf1WQXB+bGvU3KtO\nOBq/Tg2rrFBBAvaGoAESbIQu1UlGD+0iXcevxo1+9lcWgm41doYY2NSoubk6PH7NY8qLllmh\nggTsD80sdLiY2mFSNs8s1K8xF6x5xja6Mmxo1PqjXsn8Wvo1VI9eHxKwP1xRSuVoWZ/zOEH/\n2fClV17XEGk0doYY2M4oj3fHJvI7kbE+vRJ+DytgNSjSXalgiIMeeti9dvud6hdkhWsydoYY\n2MyolaVuCPzOCbgifuNt+Htch8pc16YVq7zqTABHS4YmOQh7THlc0y/XU/7odhZC9fwGLFWN\n5He6V3vdtfALkyiXoFt/5cr3DGq2/IB4CwIGkQ01xwC7cCcW/edaCIapd8VvyFLzOH4XUA+/\nLgF/PJz7O+v+jW9Y+ry65V+SbUDA0NEyZvjloTX9QvlyXlMPDVPvid/AR0Ui+HWotkJ+ZZKX\nYer6a0c78OOxKy9Ja0LAPUZ2xLhIfM7TaeuXSf1KhsFVw2D8Miz6YZwc35R74zfmUa8wft3q\nrY5flQS6TOf1xbKGjcsGtyJg2Phh8WeoX0W2dtyUvyaO8bUWgzYZJyc05b74jdRvAL9LqIlf\nf8kKYnMAACAASURBVJ5M1N5DR0D6VWLLmDPw168Zp+Sys1aemwp9+toET0YaA3ugan4Xpn8X\nEMSvW7QV8hvb+r8x0mt8pLTSMVIwGPzL5FAnUL98dNJAHpokxaw/w2zmb37UzG/cmxKC+HVJ\nt05+o1v/Aqp1tt/XglNEWcAYh/KN3N30CejXCEUDN41pJ0tk25SAK+Y38k0nAfwuarkyfo0k\nTzVQ8kj5dx/mCbvro9J5wkBAgsWXGYJPbE6/0r/ispvenmCIXfAb+6Yif37n9VshvzDJwzAV\nDY0I2BjKjF+Z08ea6lcucR/7ZMWf8NKqGQPvgt/4N4358utSbr38wiTIP7nhKqJmqIigaNUZ\nOPVrPBQ+pJZ5bR2lBNgDvwlvCvTk16nfevmFSYI65u+NscvLJ2ErAuacwd6TSYKc+mVSv0P4\nygqHcHCRgL/yIAtpaOPktKbcAb9pb/r04XdGwdXyC5Pc2WKwwsC3G6y9CiN8i2gD0jfiUxfL\nqV/V9+rvtnOEZ1g82ucX6U29C/zOCbhafo0k14v3687uvTv2fXbDMtlqCEYCDDXO6vcEpo/M\nuUE5WAq85/kYlpa8cX7jpn8dmOfXodvc/Bpv9Iqoi6qTafoauvGUT3f+VEMwFmATGDxO9Hs6\nsUnDjSl4oCflaVl8yvb5RdTvHL/rQOfXqFNBActTvpeLi+CwzKqCmiSUW7oic9O/dm2z1bmc\ngOvjF/P2O8Ovp4Qx+bV/0DgYsZac1XjqfKmkh0YCDHGY5C7pV10KnFcp4FBUxy/aD9XM87ss\n2yz82j9oHJ5DrCVPJn8c68MuNRCMBRBWHPiy5Tt7/9WeGVu54lONK4Ha+MXW75TfVQnj8zup\nVKKAgbWXxdWvPe6qEq+VDqk9ASt2QZs4pn+XiM5nXELqdvnF+6G4OX4DhYxhyrRSeAJmyw+g\n9Hhf5bfPbXOC8cDUP9VXe+jX7M6z1RlNwI3xi/hDjzP8rqk2A7+OSqW60Leun7l/deyPX5ef\nQIkuom4Avsw5wjX9mmRzOU0xXwxfPmMuYULdeLv84vnPc/yu6TcDv65KJQr4Lp4BfbML/3q8\ntSGiiKqhqOLgH9Svlu+i/6zds5liuOY40MKk+jXKL+r0kZvfMCWj8Ov+QYmIStklyy94w7h2\nBMxhsJExsb2mX3AZcDm2kqFOdzH2Z4CFKWiTXzT9zvPrIVpsft11ShRwp3roLhPBrOZJYb1I\njolIo2DN//7LxwtEz647ixGFzZ6xZGIKmuQXVb9ufr1vxEMuKPzO/SBMTK00+rcw8GGMdOf/\nVt4lGlNElNdYDopLQZLc9lq+AXrragXcIr+Y4asZfkOAxO9cpVKDWPItDMPim+VHzxx902oR\nUV5jOTCmemau2PLTLwcNUa+AG+QXUb9z/PrqFpHflR90CquVsfXqX8Jw7btp9lhO9wwneKnW\nNQCMiXS40nP5hrgoBLtMelsMZM7M+m8wBm6PX0z9uvldE3EGfhcqlSzgALw7Xx+sHQEb1A6f\n9tODs/6zNdGg98qsNd+bRaGDUAO/qPp18usNPH4Xf5EtolbhSQTevjOJzQiYA3bHDnuq33WK\nOeBUdcRMf3LNbnBLlGy67fnF1a+L3xgJJ/K7WKkEAVvehE/SJ3h1sE8RVY+BJTPaS+L6nqtD\nWTabJsUctB30ssx9KVYm1K81fvHCzz2c/M6zmY3f5UqVFXBoERVHoXWPyeRAh5k+89r9V7DO\nVQXHzESGTQrYO2tYBhq/qPp18xsKBH5X6lTUhY4pAvnSQYO8d4gOlYN7rr9+meFIjQFLmRnn\n4FCqnZUhC7/I91/xafAbI2Gexu9anaoXcIVgalQjCRkjisbv/a7oV8QgVX8s+nvAsTlGijc2\nNYMcyGEU4tNHc/zqz1U5Y/G7WqlkAffTDJxfP+H5eBdRF0yvTzpaPTdQs173X+2jgXiH7PFh\nlDLF3MT0rfCL+PQCn+HXS7mo/HpUKlXAF2Foh8pw1QIGn+LvyEygftWQiKu+XeTE1XWD0A6J\nWbTCL+7Tv05+A4HAr0+lEgX8ZJdvb4Z+GwMKKhawDuGpHnXYdXLp14N15VlxySxDrX9aVq3w\ni63fKb/aGTZFmpFfr0olCrhj3xyR4voFbDI0o98gyM4eRDgQ7Y1EI/wiP73Pp/zCKeE4BPPr\nV6lEAev6VU0wInT0UNGiQlaCY7d+XfzLiAgHq2VZTU3ZBL+400fij8mvunXOs4nNr2elEgV8\nFj30G+9hb7uI2qA7Ui4nBpjUL4djYZtR7mQVXCly7Sy+udFogV/c6SMnv6oLkz61S7Go/HrX\nCWcM/EL+EayqBSyJ5IpgufpKx6Lne2O9DScXBdEye0xzk1I3wC+yfl38KmIctObh179OiQLm\nV2Eh2qOi0yKCE2NrwFEA04FEl/+8RDRgW8cn9Zaj7kkVSmyL6vnF1i+f8qv2KGXm5jegTokC\n/hOPm/0Lz8a3iNCkOW5jjmJk16z0O1K7pl/IqPwyZAeuFdQKpbVE9fziy1cUo/kVjEjpLssX\nhd+QSqUGsboH7hT/tIiopPkFrEgCw1/hTC8ybDJtfpc5Y1YoMYhVOb+59MshKVLA46cnEvgN\nqlSigG8/sy7//H+CMqKIuJTZFKx6UMFJkH4dsUsOvztMT61QWkNUzm8G/dr8MjAc5msCRuE3\nrFKpY2D+r1+rc3uFZ+NfREzKTAI2/R2mw1c87P4LvDHgquki8CqU2hA185tj+MsNfvVfeRf2\nozeB38BKJQuY88/jzFiH9tJvVxHhKXMJ2CxERq9GWpR+nSybOzkcY4kP5xhpawHXy2+G8JXN\nr/4rHGjunEDC5De0UggC5vx7s81KRK1jYKu1VfQKxK/WumhFpmJUkSqPo1YIoyWq5Deffq2/\nkl/xf15+g+uULuB330Gzy8obz5KKCEpqOEHIsEdg2mNWW/5hDqYJHg3XQUvECiW3RKX85tDv\nVMBmKIvxzPxG1ClRwK97x9j5jjxESrvqpm2EBmsENup3pEgPho3u2EWqdKVU32u4Z8gVSmuL\navk9hd+rfLK2/nIZw5IS5Fn5jemTEgX8K/vq9xqk6CKqgo5KwuEvU940lwQ7ZhHgVKBysCbk\nYlc9Lb9a+c1y/+WQX7VjJEo7uvn4japT6h24Hx39emjkiYZ6BSz9HT17NOxRa7EAfW6YYRAw\nu6iHTLgWJ6WulN+c07/meEVJUh3Oxm9cpdLHwH+9l/UjOaZ0zyLqgeyhpX4FQWrTGaTkkE85\nr2AOiAT1vDYB8yr5zaZfxx2YwT85+Y2sVLqAf/irKEqZFSpWCATLbf1Chic9tTyHG3Tq7lle\nL2gtgJFRXfzmcp+5ewwM/2TkN7ZSCAL+9mHKcy1RyqyQsUJ1DQ1cmbdjs2+eEqwnGbgmWAy1\n5BXA9SWDY3IKKuM3o36nAtaaHLd4Ln6jK5Us4GGlzv0vsnivIuqB5S+PPaqKP8tYpXaoHM6W\nugrMHppzMNICnzgmJ6A2fjPKdypg6Q7rWFQmfuNrlSjgYa0s9iRDTQK2nB1m6ZczuH7D6p65\nxe+YH5exTMioxay+ZtYMWq9A0Nk2quMXW79TfjmHOjakmo1fUKkC/MIk9T2tggpHUNIY/kr/\nWdNpk2ozLJwsM0rJjIJUZj4GrVchotogdWX8osuX2/xywANwh3PyCzqlIvzCJMiulauILTF1\ndozbL1frN0TvK/8YoyPNOZfkwC5b/lPRTzvwuWyQXxViURm/2O6zozmVsAQxIC6diV/oVBTh\nN/6S+HuML3i4ro2pKhEws/5a7rPULxwWwR7aEeJgin/DkzJ8qwUBTwzyr0MJ5OY3j37d5ij9\ncp6X36l+c/NrJHmeQfezjO8ZVHX5FS2VCtjS77B+g5+kb6XGPUvQgyOQRoc3uZnZmkEBdYhE\nTfziPzxo/bUPaQHn49cY1JcX8AOYvoY76/6Ny/I+r275l2TrFLCt334HE8uxzL4ZRiqNPbo3\ntzpipsMazNhaMiikDnGoiN8CDx/ZR6YCRufXrFV5AYe8rbADPx77Zp1vEVvCGJJo+YpBzXD/\nZVzEKARtS3207neVJKAymKvUBYMCqhCLevjNMn1kNCegQt4iGed5+bVrVYRfmCQkYMb8E1Yj\nYOX/jMudmWpgxsT91+h8V6A9Knu/zNYudcEg3yrE1NssL+ZcbH7zTP8acSf1VfCrPOFs/E7X\npBThFya5M/9l7g3egeFQRT58xHX8ip8MT0oyN3WvxKYiWEcyOWd8QrB7AGwa5F2B0BobqIXf\nfIsnZXOCex8Dmzn5da4pK8CvkeR68Z5p+I2RXuOkYitjYA09+Wv4z3qI5Oh4Z3iGvTWDV0u4\nA+VDdmJT1sFvztVXIwAFkkq1OxO/605FJn61J2BZvIYLOPu82LPXJmBr8eSoX6YeB7afRlmA\nOkOdxb18Zwf8zo5vynr4za9feHfkoGkz8ruq32z8RguY/92HecLu+mhjHlgChq/UgFjErziY\n/LPaxNVFq33Kz+Kwun7tKM4Fn2tnxaAafgvo19QYV7Tk49fj/msYtmZ4AAqoqy4Bq5dtqB5Z\nzf+qIY4MWbJZGN3y8B0EOiNgj6xWTqsKQUaV0C83+VW+czZ+ffWbg9+DCRi8bEMOjfT8r0mw\n7V656J6GKAOjFhIHEXDWp48AIL+SlHz8etSqjICVcd1i1FGgtaWUPYxnf0cXaPSfdT8rj0pX\naxFgGCWDmnEoIuDN+S2lXw74hULNw69PrcoK+ONxHba3lJIbw1/JIbz/Si8JjJSYu2cGC3Q0\nz2GDXgu5x8BD6o35LahfrVEYP8rBr+crCbKPgV+G/es/AN3eUkpLvyOLMn6l4h12jzvZnNIM\n+voEAWsjls6Kzb4GfkvKV/LLDQFn4Nf3lSLZo9Ac9rjn9enC9hZyQPmKPngc/wLCQLRRbM9Q\nauyD/pqzaC9h+5yU0JSb81tcv1xrNhu/8G2maybl4dc5BvZJt5zQaoEKYOmXq/iVyZ48YYZM\nsR+ewMfLQiRzlOzX+XoBaQwcei4Gv4X1K11lW53yBBx+YUQUye6UJNeAt422dge29Sv8ZyZ+\njUF1teoMsTnhdvJXnsbnOlm/4Y8X0jLZkt/S+uVqzJuRX/H2Q64/EQxPSBLSizS2lNIMX41/\n9P1XDoMUVdpt8gK4UqaA10xqU5S7AyPzWzZ8Nf7h3NQjOr9Qv5vxC5OcAxa7N7WUcnL7FTv5\nSZMDPSbYCfuxq7rqKbTXxsM0tJBXHDbjt2z4Wf4Fg9cs/OoF9aq8DfiFSb4Bi91bWko5p192\n0gMlrrpmrlwlud/2s0AABPTLbIY8ZvzdUsBb8buJfrlkJBe/5usgNuPXdKFBJfCwtYCt2V+1\nc1g/KRmSn5os2ANbJHO91xhXuWEym9QaqS70JvwWnv01t7PxC2ck4dHS/O5fwMbDC5IQoV8r\n4OjskG0ownUkROQ+Y4DpWx1OwAX1yww+xj2Z+J19e+yWAs6ETQV8moavmNav7KEVd+Mu6XY5\nyeahd2BuxD82FHAmLBtVPnxlCDgTv9bbrzbkd+cCNoe/whIm41dmqzP1T7pbinbxaftZTPTR\nOo8FbD8GzoRFozYIP3M5wM3Gr/vtG/ZZMTWIT8JglXbjQhv6NeNXHA6YNMFMDpCcjpZrp+1D\nzWGjKCUofQN+i0//ir+isfPw614+uRG/sQKe1BLVKiSYLW3Fr8AeTbBk3QxtcDPsMa25nyTS\nhVNOwDj8ltUvFPBkDx6/s8ufN+E3tshnAwK29CssORlRLXlI/1EjJdn9KqoB71oQheoCLSwB\nDH6LTh8B2Zq3JY7Mb4lfRM2bZMS7W37IDKGINDhuv/0HeCOH6o6ZeY7Fqx4Oad8L7i+IgoWl\n81v46UH4kZHfIr+ImpDk78K6u9eCnffyArv5Ispgcvsd/4IXynIxWgL2QZfJuAHJGAfsoZka\nHxVDemHl+N1m+MuYsYnOb95eKUXA7x+zzx9rPTovhp9gvTuyVelw61dFtZi52wmzJ5YMc+1p\nyY1ySCirNL+b6Nfagc9vZq8iQcB/g8X3S/fm34tv35vNqmSs6NcnZAh7aMZdcUs5biqH+LJK\n87utfnPxm3tUkCDggdQ7Y/0PuH+XHx8rYFUirDjDVL/cI2SoBkpqnnASmxUBzXKIL6swvxvr\nl+fhN/uoPkHAYuzAwAYWigvYR7/ymO5p7SNcu1KKZ+F3ccUsC6nd6jW1nkN8yqL8bq9feQyR\n39Xw8yb87lDA/vqV7LqCUaJDtjwtxhTJMpjpa5aHV7eeR3zKkvzWol9Uftf0uxG/+xNwyP0X\n9LKzBBu9uKRXtlPI/dc0JwpNCLj49O/CMTx+1++/q+asgwQ8aWfFwaJ+lccEai4CG2DGELRW\nOF3M+huFFgRcVL6LfSgmv576Lc+vFrA9eEdDSQHbt1/l1jj0OyHYdIIM90o6YjplqMO0uYAL\n8Vt+9eRsbRD5XV9+RQLGgEu/gwEu/U4Jlmer5NZ0v8FSYCsdRMBbrH72FrB1tj+/HssntxZw\nRpQTsPv+y+T6qwn03IH4LllkMinsoscU3HX9+0hiozFSfkCjtrj/zs8EYPF7OtXL754EbA9/\nlSc0o18QpZT/BKt82j3LQSSbytXPnd42Cp0RWxml+Z0VMA6/J/XbO4vmbMPvfgQ8CV8pr2lO\nv1zSKfpoJneB9Lbbabha+lLgPtVMd11JwEax67FEDH5Pp5r53Y2AHfrlq/rl2jqmCVbbDHhU\npkpVbwuujd00ZSC2MQrQsbiaJplf/eMLuBVYMjZvkgqLmNHv0PC9fj0JBvdTuO24zWpCwRwF\nRj2WQQI2C5Vuaz5+xcVTL791WhWMyRor+YcZz//Ko+Z0oDxVDoHktuyZpVsFwhvT03bTlMHY\nwijNr+vmiMfvScavquW3TqsCMV0jqb449GtFG/QwynTI9MNlgmiQzDybgVR5QQK2imTm5riF\nx+948VTN7x4EPK9f8f4N5jiqCWb2hzHHYHte+kSRjbgQkgMYPiABWyW6Qr94/A7+c+X87kDA\n09WTxjG75Zn1lwMKYSxSMzyZcoAZgO47P0jA3NbrpOnx+NUvb6mY3/YFPH/7da+gmRI8bFl+\nl+6VdahDfMydUwAk4PXi0PgFv91dMb+tC3jBfV58AfeEYHOfHiaBgRNTSwfkWYb3VQAk4PXS\nsPgdxl7189u4gH31O41Lgs4YeE3MPKj/ah+LFyYVggQ8Uxg+v2UfirTsyZukoiI89Wu6QcaW\n9p1gD+wmWPlYRd0qiMML2F0WPr+O2EkJHE3Ajid84cGZH4A0+mt4Q2Xy//keejzBzq8Yji7g\nufuvdTCZX+PdwwVxLAG7ntB3HnWPiriOQMrvIuTIHZfEOEbSPtY2XvTBBbyoX0R+R/22wW+7\nAvbWr5tgHXwUzpMKYDB9FJwNPSxxrDzFxxbwXEHY/Mq1A03w26yA/fU7Q7D4ZFx7UMYRiz1j\n4pCrZGVxaAHPloPMr/Kfm+C3VQEH6Bc6TDaLgGHO9Yy+nbWc6JddNdBwURxZwAvFoPIrll+1\nwm+jAna84Nk46libZUaPFcFmlHLqXClvWd1/FdW8BYLzo4hRi94sJr99+LklfpsU8PLt17V+\ngwE6TIJl18s5/G4QDE8U91+dpCgOK+C1MtD4FU8vtMNviwIO1695lsGw9cnG/8GkInC+FLdb\nzSQdVcB+RSDwK5ZfNcRvNa3vj1j9WgTrIROzCNYRSPlND5J0WFMnLIeDCtizhGR+h+UbjfHb\nnoBd73c2Ds+ugWPmXwY4Vrv14Ef6ZAz+B66F8pMMRxWwbwGp/OrlVw3x25qAne9nXz5unzln\nDyRcc6wZNvdsgkMK2D//NH71qx9a4rcxASfpd9Ev0hP5ZietfCwdAmEtEZwfmY0KyD6J39F/\nbo7ftgScpl++5BcxdRgELUEqFYHeYnrBMLI25DUqLPd4fkf/uT1+4y39e1yHi/p6/8tVhI1k\n/c4DDpSY/Q04X8teWm6ULHcDfjNmvsavil+hlhqKggL+nvXcN7tkKcLG0us1Zk/whUHw5ENN\nAbMlLy0/ypW7Ab9Z817hV7jPDfIba+qddf/ew7fPq2P3HEVYcC7PWDnBH2YE05gnZGAP4wte\nWn6UK7k8v66c8bJe5lfcflvkN9bYjr3V9zfrchRhYs19TtOv7JrZlD5zwASObsB0uQKL85s7\n4yV+1ctjG+Q31sbJalL8Igzk1i9YNGvVZkbAm/ha5YorzW/2fOf5FU8fNcpvG3fg1eFvsn65\n9KaYnfOcgB1GZMeR7sDo2c7wO145zfKbMAZ+fYZvBcZIZfQr8mTK3RKfkEqL39IMFx0Dl+O3\nSK4z/Morp1V+ow28gCjl+ZulCIki+jUJZqaG9addfvUEx6Igv2UydfMrl082y2+8gX/3YZ6w\nuz4yzxOW0a9FsC5FRT6MoEYzBEejGL+F8nTyK5ZfNcxvAQMTiyikX8mwwS9Tz6vMmFD/GCk/\n8ogtS6Y2v+OPhzbNb64IBERKRos/rrBwSgRAlBIImGvC7bN5C1HKPEDj1507eo5DrhN+1Y//\nNsxvtIWfG+senD/PrFsMcaTR4bz9Jup3/pKzrkk548DcdWhjnjAWZfh15peYoTe/K7/92wa/\nsTZ+u74dno+hObIttVu//Qbrd7FntaYKN3u52RzKmVGG3wzZ+fO74cvr5lBQwPd+auHesduX\nf++5phky6Hd5bKOHSSKqob6GFZIL5cwowm+O3Lz5lU8PNs5vrOWdiNAOEwx5Jvqd0kTR79zA\nDQYqjchkJfwWtKMAv1ky8+ZXDH+b5zfWdLnEH2wgF5FFv+DuOnMQ/uVLJ2+Bcnbk5zdPXp78\nqguneX5T78D95zdHD51Hv5pDTwFvEsuYRfk7cP+Zhd9MWfnxa/zwTtv8po6B71/xHbmITPpV\ntxa3WQz8XyPKj4Fz8ZstJx9+sSYe0VFQwHmjlO4Gnug3PGPgMzkJdj+wUg3KmbVFFBolIw9+\nxfRvhSgo4KzzhBn1qyMcbrMyLU1AQkHDys8DI+Wzym+199+yAs5YRC79amahSBf0Wp2YqzJG\nAkt46Tl48Qt/OHoH/FYoYPfwN1m/2reC65sXopDCma6J4opM0cDzfNNzWOfXlG/7/FYnYK/b\nb9T9V32KCXxrryMBU5+VoCJTNFBCT1hZrPFr/O77HvitTcCZ9QsJY8ZeRwJmJdge9VgCgCY+\nlCyW+Z3ot3l+KxNwNv2SgLMBYfCKaMUiv/DC2Qm/dQk4n35JwNmQPnrFtGKJX+PC2Qm/VQnY\nLUwU/RpjJGs7YAy8ccijnksNINUopEqt8mtPH+2D34oEPCNMJP2CeCSMTAZFoTdfOLtHAWPV\naY3fyfTvPvitR8CZ9cth92qyNmugPU+4cLsugx0KGLFKi/y6Xw3RPr/VCHhGl4j6jYFcG2A8\nNLodw/sTcKkazSy/ap/fWgRcpX7BY/3DMtoZowoatFnJC0gxalv97oHfOgQ85z5vff+FVjRK\ncH7EG1VsuDl3/4VWNMpvFQL2HP4i6Hc5ymgdZXKv3Bx//pmt5JIR+xJwhtq4mdFP7++P3xoE\n7Ok+p+t3Oco4OTolWLhcS7lkxa4EjF+ZGWbkj6fskt8KBFxMvytRxsnR2R56KZes2JOAc9x/\nnRmr+699dBf8bi7gOVlm0++cQY6jzjHSlgzvSMDZ9GtlfTL1uzt+txZwQf1GCNgVpWyM4PyI\nMipHTVz86vDVTvndWMBzqsyhX/2M6GKRZpxjOk/YGMH5EWNUloo4+IVP77uKbp/fbQVcUr9M\n/VbZYqFLDcI8z8uGvQg4y/3Xwa8xfbRPfrcU8Kwq89x/AcfuE1ajjy1GKfMj3Kg8998pv+b0\n7z753VDAs8NfPP3CVyPZO5ZPXzyjpXnC/Ag1CrHxlvl1Pb7glWFL/G4nYO/bL8LzRyDbKjXg\niyqNDzQKU758id+K3z45h5YEnF+/5mCGBJwLYUYhVmGZ3wb125KAi+nXYrhKCXijSuuDjELX\n7xy/Dcq3IQH7h68SiLAJ3i42gYYqrQ8xCrMCi/w2qd9mBDzbuJj6nTpVVb0COApV2h9gFKr9\nS/y2qd9WBFxGv/twmk1UWRd/o5DNn+W3xeHvgDYEXEq/u3CaTVRZF2+jsK2f47dZ/TYh4PnG\nxdYv34PTbKLK2vgalcF4J7/t6rcFAc+3bQb92iXAha9NokrDPY3KbzuTb5+sspl8UL+A52+/\n2fVrPXrSJKo028uo/C0uiB3uvwfid7uVWGvnoN9/x8+mI1tVmh3JLzbGIoT/XGVDraNZARfT\n72TlXVuo0uo4fvNYMXl8vy20KuAC+iUBZ0MUv3ms0OPfKltqFY0KuIR+ScDZEMNvHit+V83h\n+K1AwK7jGSYC5BgY/Ppzc6jS7Bh+M4DxMXx1MH63F3Ah/cIodLMTSVWaHcFvDjCp32Pxu7mA\nS+mXc/nOUNYsvy0KuJjJv+HvEfmNr+rf4zr4K9f7X0IRZfTL7KFRowyXNDsfv/hgjJnTR8fh\nN7am3zPTuEQXUUS/YO3G8QiOREZ+0cGYCj8fj9/Ymt5Z9+89fPu8OnaPK8L9QzaRFq0ZcFCC\nI5GPX3wM+j0qv7E17dhbfX+zLqoI54Fc+jUU3Ci/Be3Oxi8+mFi+cUx+Y6tqdK7LPe3cwUL6\ntQSs3ekGUc7uXPxmAIPLJw/H73Z34FL6td2qdmOUbd2BC4afT7DAg/GbMAZ+fYZvkWOkYvpt\n3a2CKDoGzsBvDlj6bRolp5EuIEp5/gYXUVC/jbtVEAXrkIPfHJDLJ4/Kb8I88H2YJ+yuj4h5\nwpL65W27VRBF54HR+c0BMfw9Lr/brMQqrN/doMrL1JffHGj37TlO1CNgBuFX6L6oyINqBBzB\nbw7sTL9FBfy996HJx5mxy7+wItzuzjZUNOZ7FTQWm98cWNfv/vmNrd+n+zXNt4tYaucucRP9\nNhf9KGcqNr85sKrfI/AbW7sbu35/H7fPj+tbyDRDRfptb36pnKnI/OaAx/0XfDaBoiux1p5j\n+wAACTpJREFUvuLj520FTPTXp9+WGC65EguT3xzw1e/O+U1aStkxsOFVRE36PQbBsSVh8psD\nHvGrQ/Ab70K/OX+M6+2+y4Mk5vwKsFUo8RAERwKT3xzwiT8fgt/Yyr1Zd3/za/dj+HVmL78i\n6tLvMcZIkUDkNwf85o+OwG907V6dngh8+BVRm34PEaWMBRq/OeB5zRyB34Ta/bsNb224Pj5+\nRVSnX36EecJ44PCbAQHLN/bPb7mllLn02xhHSaiypov8ImRv87u35VcAFQt4RmXJXDTnJSWh\nynou8Zue+YTfHeu3YgHPFINw/wWf+0eV9VzgFydzI/s967deAWfWb6VXNj6qrOY8vyh5G3/3\nrd9qBZxNvyTgCjDLL0re5t9967dWAefTLwm4AuSMQFj87ly/tQrYvRuHCxoDb46sRhn87l2/\nTQkYiQuKQm+OvAIG/O5evy0JGI8LmgfeGJmNUvzuX78NCXj/XGTBEQUscYRLphkBH4GMHDiu\ngA9w++XtCPgQZOTAYQV8DP22IuBjkJEDRxXwQfTbiIAPQkYOHFTAR9FvGwI+Chk5cEwBH0a/\nTQj4MGTkwCEFfBz9tiDg45CRA0cU8IH024CAD0RGDhxQwEfSb60CJqAhP1vhyFrh0ylr9rUh\novXxCS1X4t6zqlKva6iyJevMCiMnEnDFWZGA950VCXjnWZGA950VCXjnWZGA950VCXjnWZGA\n950VCXjnWZGA950VCXjnWZGA950VCXjnWZGA950VCXjnWZGA950VCXjnWZGA950VCXjnWZGA\n950VCXjnWZGA951VmwImEAhoIAETCA2DBEwgNAwSMIHQMEjABELDIAETCA2DBEwgNAwSMIHQ\nMEjABELDIAETCA2DBEwgNAwSMIHQMEjABELDIAETCA2DBEwgNAwSMIHQMMoI+CmKgT/hdO9Y\nd/+GZnRWiUD6xKxSrfreGLu9OYJVIKfkpioJ4jc4Jyx+iwj4LQx9A6svw7dzWEb3IVH3NdMn\nZpVsVTckelvpY7LSOSUbVRLEb3BOaPyWEPC7UwRf5b4/1r37A39BGbHbt+/ub0b61KxSrbr3\nmdyHTBKtAjmlGlUSxG94Tmj8FhDwk10EwU/2kDvv7PX7/Kd3+OA6ZtPnBtKnZpVqVce+IqdU\nq0BOqUYVBPEbkRMavwUEzO5cEfyUO6/sw41+KCRDZqRPzQrHKtZxJKuGnJCaqgSI34ic0Pgt\nIOA3lwRf2ev2G64P5eoeMhRfdjHSp2aFYtV9IATDqjEnnKYqAuI3Iic0fstcFYrgAReeZPWz\ndzlwCB6yQrDqH2M2F5FZyZxwmqoUiN/QnND4LSpgxv79esWhC4q3+tNdORLBMqtkq57XbhjC\npFulc0JoqmIgfsNzQuK3qIBHfPuIebTV3+4CEiYRLLLCsIrzm8VFfFY3PTxKNaoUiN/QnHCM\n4psIeNjqYq2+jPNlIH1qVhhW9Vx0OFaNOeEYVQrEb2hOOEbx7QQ8ht4+oaG3z/nyGb6A9KlZ\npVvlSp+QFWitVKMKgfgNzQnJqMICHufBBlMfw+TXaxzRe+PFpFME0qdmlWqVTH9OtgrklNxU\nRUH8huaExm9RAd97I7/DxHXU8pOPIiV5pQ7IKtWqYX3N99qPbDBW6ow5pRpVFsRvaE5o/BYV\n8HdcCzp0NWcVRvfHjekVpCB9YlapVokVrrYpMVnpnJKNKgriNzQnNH7LjoG/946dn+prF+g1\nMM0KTI+RVYJVw9Mk0/RRWZk5JRlVEsRvVE4Y/FYc2iQQCGsgARMIDYMETCA0DBIwgdAwSMAE\nQsMgARMIDYMETCA0DBIwgdAwSMAEQsMgARMIDYMETCA0DBIwgdAwSMAEQsMgARMIDYMETCA0\nDBIwgdAwSMAEQsMgARMIDYMETCA0DBIwgdAwSMAEQsMgARMIDYMETCA0DBIwgdAwSMAEQsMg\nARMIDYMETCA0DBIwgdAwSMAEQsMgARMIDYMETCA0DBIwgdAwSMAEQsNoX8Df55Uxdn0unPJy\n7mAedZ+c45OIgAnidxGt2TvBq2Mjus/cKWe7kuOOYxDcOojfZbRmr40XY7e/39+/K+vmzklg\nqX2CGwfxu4LW7LXRMek/3dicl3VoghsH8buC1uy18I9d5dfPve+p+fPMupHp14Wxy6vnhA20\nvH5jqe7O1Y6Rq9/55+F8xj6/Xv7B4ZkGn/eO3bmR05edhyPDX1kcARHE7xoaF/CV/dk7elx+\n357j0Okp+XyM23eT4Is6/8dZ//UBz4QED2dezZxE8f9+qVRxBEQQv2toXMC2x/Nily//Xnq/\nq2PvvunP8iTG/vXbjOsd/Xb35u+uP8T6pM/xfOvMHvJMBo+/2K0/dmMfUBwBD8TvGvYhYBGo\n7LvML+9dnmu/72WepL9rgq/DSa++i2ZDb6vPtQgeO+OXdfw8lDdeFfW5V+2D+F3DzgTMmPp6\n/zlE7zc4iX9ej4tFsEoPd07ONEoyjj97l+yv/9DFEfBA/K6hcQHDMZJFMH90YvZQcHNRR9YI\nts/k3CRYH//2kxsP9oHFEfBA/K6hcQH/GwcpAwBhAq/7WY+Rbuz8fH18CJ6cyblxJjj+65hf\n/Hw2iiPggfhdQ+MCBvOEXz3mAbD6XptgOUa6Tli3CR7P/LOO8ze7vNnDKI6ACOJ3BdUZFIjX\nb2gyrNT5jVE6EUv8jVyuffzhn4pS9o5PH8R4y5HPuMOKUnIuCQZn6oJAlFId74vphtx0cQRE\nEL8raF3A/E+ulR2drXH40g9V/o17//q277m/M2sHGO6M84RcfIIzQY87zEDerOPDFTaQ+g/s\nI+CB+F1G8wL+Ne31x/HlIcILzx97t+H7sHRm6L3PwzraW785eFPjjpG8Z6dW6qhPfSZ0mR5q\npY4+Pnh2o4eliiPggvhdwg4EvC1eo79G2Clq55cEnIhLfavrCIionV8ScBLE+IqwU9TPLwk4\nCZ1+WoawQ9TPLwmYQGgYJGACoWGQgAmEhkECJhAaBgmYQGgYJGACoWGQgAmEhkECJhAaBgmY\nQGgYJGACoWGQgAmEhkECJhAaBgmYQGgYJGACoWGQgAmEhkECJhAaBgmYQGgYJGACoWGQgAmE\nhkECJhAaBgmYQGgYJGACoWGQgAmEhkECJhAaBgmYQGgYJGACoWGQgAmEhvEfna+LbXf6/+UA\nAAAASUVORK5CYII=",
      "text/plain": [
       "Plot with title \"95% Prediction intervals\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Set the graphical space so that two plots are shown side-by-side in one row\n",
    "par(mfrow=c(1,2))\n",
    "options(repr.plot.height=5, repr.plot.width=8)\n",
    "\n",
    "#Confidence intervals for predicted values\n",
    "plot(data$Gestational.Days, data$Birth.Weight, xlab=\"Gestational days\", ylab=\"Birthweight (oz)\", main=\"95% Confidence intervals\")\n",
    "abline(model1)\n",
    "\n",
    "conf_interval<-predict(model1, newdata=data, interval=\"confidence\", level=0.95)\n",
    "lines(data$Gestational.Days, conf_interval[,2], col=\"blue\")\n",
    "lines(data$Gestational.Days, conf_interval[,3], col=\"blue\")\n",
    "\n",
    "#Reference ranges\n",
    "plot(data$Gestational.Days, data$Birth.Weight, xlab=\"Gestational days\", ylab=\"Birthweight (oz)\", main=\"95% Prediction intervals\")\n",
    "abline(model1)\n",
    "\n",
    "conf_interval<-predict(model1, newdata=data, interval=\"prediction\", level=0.95)\n",
    "lines(data$Gestational.Days, conf_interval[,2], col=\"blue\")\n",
    "lines(data$Gestational.Days, conf_interval[,3], col=\"blue\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, we see that the prediction interval is much wider. Loosely speaking, the plot on the left shows a range of uncertainty about where the *average* line is. The plot on the right shows a range of uncertainty about where individual measurements will lie."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}